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Abstract 

We review a resistor network approach to the numerical solution of the inverse problem of electrical 
impedance tomography (EIT). The networks arise in the context of finite volume discretizations of the 
f~^ elliptic equation for the electric potential, on sparse and adaptively refined grids that we call optimal. 

Cn The name refers to the fact that the grids give spectrally accurate approximations of the Dirichlet to 

"■^ Neumann map, the data in EIT. The fundamental feature of the optimal grids in inversion is that they 

1—5 connect the discrete inverse problem for resistor networks to the contiimum EIT problem. 

^ 1 Introduction 






X 



We consider the inverse problem of electrical impedance tomography (EIT) in two dimensions [TT] . It seeks 
the scalar valued positive and bounded conductivity cr(x), the coefficient in the elliptic partial differential 
equation for the potential u £ H^{il), 



, V • [ct(x)Vw(x)] = 0, xen. (1.1) 

> 

^\i The domain f2 is a bounded and simply connected set in M^ with smooth boundary B. Because all such 

fy\ domains are conformally equivalent by the Riemann mapping theorem, we assume throughout that fi is the 

"^^ unit disk, 
1^ f^ = {x= (rcos6i,rsin6'), re [0,1], 61 e [0, 27r)} . (1.2) 

, The EIT problem is to determine ct(x) from measurements of the Dirichlet to Neumann (DtN) map A^- or 

!" equivalently, the Neumann to Dirichlet map Aj.. We consider the full boundary setup, with access to the 

entire boundary, and the partial measurement setup, where the measurements are confined to an accessible 
subset Ba of B, and the remainder Bj — B\Ba of the boundary is grounded {u\bj = 0). 
^ The DtN map A„ : H^^^{B) -^ H^^^^{B) takes arbitrary boundary potentials u^ in the trace space 

H^^^{B) to normal boundary currents 

AcrUB(x) = o-(x)n(x) • Vw(x), xeB, (1.3) 
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where n(x) is the outer normal at x E B and u(x) solves (1.1 1 with Dirichlet boundary conditions 



u(x)=Mg(x), xeB. (1.4) 

Note that A^. has a null space consisting of constant potentials and thus, it is invertible only on a subset J^ 
of iJ^i/2(6), defined by 



J = Ij e H-^/^{B) such that / J(x)ds(x) = ol 



(1.5) 



Its generalized inverse is the NtD map K\ : J ^ H^/^{B), which takes boundary currents Jg £ J to 
boundary potentials 

Atjg(x)=w(x), xeB. (1.6) 



Here u solves (1.1) with Neumann boundary conditions 

cr(x)n(x) • Vu(x) = Jg(x), x£B, (1.7) 

and it is defined up to an additive constant, that can be fixed for example by setting the potential to zero 
at one boundary point, as if it were connected to the ground. 

It is known that A^. determines uniquely a in the full boundary setup [S] . See also the earlier uniqueness 
results [Ml HH] under some smoothness assumptions on a. Uniqueness holds for the partial boundary setup 
as well, at least for a G C'^+'^(J7) and e > 0, [^. The case of real-analytic or piecewise real-analytic a is 
resolved in [371 [Ml iSl iS] • 

However, the problem is exponentially unstable, as shown in [H [51 113|. Given two sufficiently regular 
conductivities ci and (72 , the best possible stability estimate is of logarithmic type 

llcTl - cr2||L°-(n) < c|log||A^i - ActJ|^i/2(b)^j^-i/2(b)| , (1.8) 

with some positive constants c and a. This means that if we have noisy measurements, we cannot expect 
the conductivity to be close to the true one uniformly in Q, unless the noise is exponentially small. 

In practice the noise plays a role and the inversion can be carried out only by imposing some regu- 
larization constraints on a. Moreover, we have finitely many measurements of the DtN map and we seek 
numerical approximations of a with finitely many degrees of freedom (parameters) . The stability of these 
approximations depends on the number of parameters and their distribution in the domain fl. 

It is shown in [2 that if a is piecewise constant, with a bounded number of unknown values, then the 



stability estimates on a are no longer of the form (1.8 1, but they become of Lipschitz type. However, it is 



not really understood how the Lipschitz constant depends on the distribution of the unknowns in f2. Surely, 
it must be easier to determine the features of the conductivity near the boundary than deep inside fl. 

Then, the question is how to parametrize the unknown conductivity in numerical inversion so that we 
can control its stability and we do not need excessive regularization with artificial penalties that introduce 
artifacts in the results. Adaptive parametrizations for EIT have been considered for example in |431 150] and 
[21 [3] • Here we review our inversion approach that is based on resistor networks that arise in finite volume 



discretizations of (1.1 ) on sparse and adaptively refined grids which we call optimal. The name refers to the 



fact that they give spectral accuracy of approximations of A^ on finite volume grids. One of their important 
features is that they are refined near the boundary, where we make the measurements, and coarse away from 
it. Thus they capture the expected loss of resolution of the numerical approximations of a. 

Optimal grids were introduced in [29l [30l STJ [7[ [6] for accurate approximations of the DtN map in forward 
problems. Having such approximations is important for example in domain decomposition approaches to 
solving second order partial differential equations and systems, because the action of a sub-domain can be 
replaced by the DtN map on its boundary [61j . In addition, accurate approximations of DtN maps allow 
truncations of the computational domain for solving hyperbolic problems. The studies in [2H EHl SD H H] 
work with spectral decompositions of the DtN map, and show that by just placing grid points optimally in 
the domain, one can obtain exponential convergence rates of approximations of the DtN map with second 
order finite difference schemes. That is to say, although the solution of the forward problem is second order 
accurate inside the computational domain, the DtN map is approximated with spectral accuracy. Problems 
with piecewise constant and anisotropic coefficients are considered in |31[ |B] • 

The optimal grids are useful in the context of numerical inversion, because they resolve the inconsistency 
that arises from the exponential ill posedncss of the problem and the second order convergence of typical 



discretization schemes applied to equation (1.1), on ad-hoc grids that are usually uniform. The forward 
problem for the approximation of the DtN map is the inverse of the EIT problem, so it should converge 
exponentially. This can be achieved by discretizing on the optimal grids. 

In this article we review the use of optimal grids in inversion, as it was developed over the last few years 
in [m [T31 [131 1371 [ini HH [51]. We present first, in section [3J the case of layered conductivity a — a{r) and 
full boundary measurements, where the DtN map has eigenf unctions e'*^^ and eigenvalues denoted by f{k'^), 
with integer k. Then, the forward problem can be stated as one of rational approximation of /(A), for A in 
the complex plane, away from the negative real axis. We explain in section [3] how to compute the optimal 
grid from such rational approximants and also how to use it in inversion. The optimal grid depends on the 
type of discrete measurements that we make of A^. (i.e., /(A)) and so does the accuracy and stability of the 
resulting approximations of a. 

The two dimensional problem a = a{r,9) is reviewed in sections H and pj The easier case of full access 
to the boundary, and discrete measurements at n equally distributed points on B is in section |4] There, the 
grids are essentially the same as in the layered case and the finite volumes discretization leads to circular 
networks with topology determined by the grids. We show how to use the discrete inverse problem theory for 
circular networks developed in [551 1131 HOI HSl [55] for the numerical solution of the EIT problem. Section [s] 
considers the more difficult, partial boundary measurement setup, where the accessible boundary consists of 
either one connected subset of B or two disjoint subsets. There, the optimal grids are truly two dimensional 
and cannot be computed directly from the layered case. 

The theoretical review of our results in [TH [T31 [131 133 [13 [111 ISl] is complemented by some numerical 
results. For brevity, all the results are in the noiseless case. We refer the reader to [TT] for an extensive study 
of noise effects on our inversion approach. 
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Figure 1: Finite volume discretization on a staggered grid. The primary grid lines are solid and the dual 
ones are dashed. The primary grid nodes are indicated with x and the dual nodes with o. The dual cell 
Cij, with vertices (dual nodes) Pi±i j±i surrounds the primary node Pij. A resistor is shown as a rectangle 
with axis along a primary line, that intersects a dual line at the point indicated with D. 

2 Resistor networks as discrete models for EIT 



Resistor networks arise naturally in the context of finite volume discretizations of the elliptic equation (1.1 ) 



on staggered grids with interlacing primary and dual lines that may be curvilinear, as explained in section 2.1 



Standard finite volume discretizations use arbitrary, usually equidistant tensor product grids. We consider 
optimal grids that are designed to obtain very accurate approximations of the measurements of the DtN 
map, the data in the inverse problem. The geometry of these grids depends on the measurement setup. We 



describe in section |2.2| the type of grids used for the full measurement case, where we have access to the 
entire boundary B. The grids for the partial boundary measurement setup are discussed later, in section [5J 



2.1 Finite volume discretization and resistor networks 



See Figure Il| for an illustration of a staggered grid. The potential u(x) in equation ( 1.1 ) is discretized at the 
primary nodes Pi.j, the intersection of the primary grid lines, and the finite volumes method balances the 
fluxes across the boundary of the dual cells Cij , 



/ V • [cr(x)Vu(x)] dx = / CT(x)n(x) • Vu(x)ds(x) = 0. 

JCi,j JdCij 

A dual cell Cij contains a primary point Pij, it has vertices (dual nodes) Pi±ij±i, and boundary 



(2.1) 



9C,j = Eij+i U E,^i , U E 



i+JJ - -hJ-i '^ ^i-jj' 



(2.2) 



the union of the dual line segments E, ,-, i = (P,_i ,+ i,P,_l1 ,+ 1) and E, , 1 ,- = (P.+ i ,_i,F,+ i ,ii). Let 
US denote by T' = {Pi,j} the set of primary nodes, and define the potential function U : V ^ M. as the finite 
volume approximation of u(x) at the points in V, 



U,,j « u(P^j), P^J e V. 



(2.3) 



The set V is the union of two disjoint sets Vx and Vb of interior and boundary nodes, respectively. Adjacent 
nodes in V are connected by edges in the set f C P x P. We denote the edges by E^j^i = {Pi,j,Pij±i) 
and£;,±ij = {Pi±i^j,Pij). 

The finite volume discretization results in a system of linear equations for the potential 

7,+ i,j(t/^+ij - (7,,,) +7,_i,^-([/,_ij - C/.j) +7z,j+i(C^^,j+i - U,^J)+l^,J-^(U^,3-l - U^,j) = 0, (2.4) 
with terms given by approximations of the fluxes 



/ cr(x)n(x) • Vu(x)ds(x) « T,j±i(C^ij±i - C^^,i), 
I cr(x)n(x) • Vu(x)ds(x) K 7»± 1 j(t/i±ij - f7»j). 



(2.5) 



Equations (2.4) are Kirchhoff's law for the interior nodes in a resistor network (r,7) with graph T = {V,£) 
and conductance function 7 : £" — )■ M+, that assigns to an edge like E^^i^ in £ a positive conductance 7^. 



1i±} 



At the boundary nodes we discretize either the Dirichlet conditions (1.4), or the Neumann conditions (1.7), 



depending on what we wish to approximate, the DtN or the NtD map. 

To write the network equations in compact (matrix) form, let us number the primary nodes in some 
fashion, starting with the interior ones and ending with the boundary ones. Then we can write V = {Pg}, 
where Pg are the numbered nodes. They correspond to points like Pij in Figure fl] Let also U^- and U^ be 
the vectors with entries given by the potential at the interior nodes and boundary nodes, respectively. The 
vector of boundary fluxes is denoted by Jg. We assume throughout that there are n boundary nodes, so 
The network equations are 



Ug,Je e 



KU = 




U: 




K 



K. 



K. 



■"^IB "-BB 



(2.6) 



where K = {Kij ) is the Kirchhoff matrix with entries 



K, 



-l{E), 


iii^ j and E ^ {Pi,Pj) 


0, 


ifi^j and (p^,Pj) ^ £, 


E ^(^)' 


if i = j. 


k: E=(pi,pk)G£ 





e£, 



(2.7) 



In (2.6) we write it in block form, with K^^; the block with row and column indices restricted to the interior 



nodes, K^^g the block with row indices restricted to the interior nodes and column indices restricted to the 
boundary nodes, and so on. Note that K is symmetric, and its rows and columns sum to zero, which is just 
the condition of conservation of currents. 

It is shown in |22j that the potential U satisfies a discrete maximum principle. Its minimum and maximum 
entries are located on the boundary. This implies that the network equations with Dirichlet boundary 
conditions 

(2.8) 



K,,U, 



-K.gUg 



n = ri 




il/2 








Wi/2 



Figure 2: Examples of grids. The primary grid lines are solid and the dual ones are dotted. Both grids 
have n = 6 primary boundary points, and index of the layers £ — 3. We have the type of grid indexed by 
"^1/2 — on the left and by mi/2 = 1 on the right. 

have a unique solution if K^-g has full rank. That is to say, K^;^; is invertible and we can eliminate XJ^ from 



(2.6) to obtain 



J, = (K,, - K,,K;^^K,,) U, = A^U,. 



(2.9) 



The matrix A^ G M"^" is the Dirichlet to Neumann map of the network. It takes the boundary potential 
Ug to the vector Jg of boundary fluxes, and is given by the Schur complement of the block K^g 



A^ — -'^BB ^BT^XX T:tS- 



(2.10) 



The DtN map is symmetric, with nontrivial null space spanned by the vector Ig e M" of all ones. The 
symmetry follows directly from the symmetry of K. Since the columns of K sum to zero, Kl — 0, where 1 
is the vector of all ones. Then, (2.9) gives Jg = = A-yle, which means that Ijg is in the null space of A-y. 



The inverse problem for a network (r,7) is to determine the conductance function 7 from the DtN 
map A^. The graph T is assumed known, and it plays a key role in the solvability of the inverse problem 
[221 [531 HOI HSl [IS] ■ More precisely, F must satisfy a certain criticality condition for the network to be 
uniquely recoverable from A^, and its topology should be adapted to the type of measurements that we 
have. We review these facts in detail in sections [SJH We also show there how to relate the continuum DtN 
map A(j to the discrete DtN map A^. The inversion algorithms in this paper use the solution of the discrete 
inverse problem for networks to determine approximately the solution cr(x) of the continuum EIT problem. 



2.2 Tensor product grids for the full boundary measurements setup 

In the full boundary measurement setup, we have access to the entire boundary B, and it is natural to 
discretize the domain (1.2) with tensor product grids that are uniform in angle, as shown in Figure^ Let 



Mj - 1) 



n 



2n {j - 1/2) 



J = l, 



,n, 



(2.11) 



be the angular locations of the primary and dual nodes. The radii of the primary and dual layers are denoted 
by Tj and f^, and we count them starting from the boundary. We can have two types of grids, so we introduce 
the parameter ?7ii/2 G {0, 1} to distinguish between them. We have 



when ?7ii/2 = 0, and 



ri = n>r2>?2> ■■■>ri>?i> rt+x > (2.12) 



ri = ri > r2 > r2 > ■ ■ ■ > ri > f^+i > r^+i > (2-13) 



for mi/2 = 1. In either case there are £ + I primary layers and £ + mi/2 dual ones, as illustrated in Figure 
[2J We explain in sections Is] and [4] how to place optimally in the interval [0, 1] the primary and dual radii, so 
that the finite volume discretization gives an accurate approximation of the DtN map Ag- . 

The graph of the network is given by the primary grid. We follow |221 123| and call it a circular network. 
It has n boundary nodes and n(2i + mi/2 — 1) edges. Each edge is associated with an unknown conductance 
that is to be determined from the discrete DtN map A-,, defined by measurements of A^r, as explained in 
sections p^ and PI Since A^ is symmetric, with columns summing to zero, it contains n{n — l)/2 measurements. 
Thus, we have the same number of unknowns as data points when 

72 — 1 

21 + 1711/2 ^ 1 = 1 "■ — odd integer. (2-14) 

This condition turns out to be necessary and sufficient for the DtN map to determine uniquely a circular 
network, as shown in [26l [23l [13] . We assume henceforth that it holds. 

3 Layered media 

In this section we assume a layered conductivity function cr(7') in fi, the unit disk, and access to the entire 
boundary B. Then, the problem is rotation invariant and can be simplified by writing the potential as a 
Fourier series in the angle 9. We begin in section |3.1| with the spectral decomposition of the continuum 
and discrete DtN maps and define their eigenvalues, which contain all the information about the layered 
conductivity. Then, we explain in section |3.2| how to construct finite volume grids that give discrete DtN 
maps with eigenvalues that are accurate, rational approximations of the eigenvalues of the continuum DtN 
map. One such approximation brings an interesting connection between a classic Sturm-Liouville inverse 
spectral problem [Ml [THl [Ml [SH jSS] and an inverse eigenvalue problem for Jacobi matrices ^ , as described 
in sections |3.2.3| and |3.3[ This connection allows us to solve the continuum inverse spectral problem with 
efficient, linear algebra tools. The resulting algorithm is the first example of resistor network inversion on 



optimal grids proposed and analyzed in ^14,, and we review its convergence study in section 3.3 



3.1 Spectral decomposition of the continuum and discrete DtN maps 

Because equation ( |1.1| is separable in layered media, we write the potential u{r,d) as a Fourier series 

u{r,e)^VB{0)+ Y. «(^fc)e''': (3-1) 



with coefficients v{r, k) satisfying tlic differential equation 



r d 
<7(r) dr 



ra{r) 



dv{r, k) 



dr 



k^v{r,k)^0, re (0,1), 



and tlie condition 



v{0,k) = 0. 



Tlie first term ub(0) in (3.1 1 is the average boundary potential 

r2Ti 



VBiO) 



1 

2^ 



u{i,9)de. 



(3.2) 



(3.3) 



(3.4) 



The boundary conditions at r = 1 are Dirichlet or Neumann, depending on which map we consider, the DtN 
or the NtD map. 

3.1.1 The DtN map 

The DtN map is determined by the potential v satisfying (3.2||3.3), with Dirichlet boundary condition 



v{l,k) = we(fc). 



(3.5) 



where VQ{k) are the Fourier coefficients of the boundary potential uts{6). The normal boundary flux has the 
Fourier series expansion 



a(l) 



du{i,e) 

dr 



= K,u,{e) = a{l) J2 



feGZ.fe^O 



dw(l,fc) 
dr 



ike 



(3.6) 



and we assume for simplicity that (t(1) = 1. Then, we deduce formally from (3.61 that e are the eigen- 
functions of the DtN map A^, with eigenvalues 



dr 



/v{l,k). 



(3.7) 



Note that /(O) = 0. 

A similar diagonalization applies to the DtN map A^ of networks arising in the finite volume discretization 



of ( 1.1 1 if the grids are equidistant in angle, as described in section 2.2 Then, the resulting network is layered 



in the sense that the conductance function is rotation invariant. We can define various quadrature rules in 



(2.5), with minor changes in the results [121 Section 2.4]. In this section we use the definitions 



7., 



'^+^-'^ z{r,+,)~z{r,) a, 
derived in appendix [A[ where hg — lixju and 

z{r) -- 



i3,q^ 



z{?j+i) - z(fj-) ^ aj_ 
ho ha' 



(3.8) 



^ dt 



ta{ty 



z{r) . /" ^-fdt. 



(3.9) 



The network equations (2.4 1 become 



1 (Uj + l,q-Uj^q Uj,q-Uj-l^q\ 2Uj^q - Uj^q+l - Uj^q- 



ai 



a,_i 



hi 



and we can write them in block form as 

1 fV.+i-V, V,-Vj-i 



a, 



"7-1 



[-d'o] U, - 0, 



where 



Uj — (Uj,!,. ■ ■ ,Uj^n) 



],nl , 



(3.10) 



(3.11) 



(3.12) 



and [— 9g] is the circulant matrix 



[-91] - 



/ 2 -1 



hi 



-12 10 



V -1 



-1 \ 


-1 2 / 



(3.13) 



the discretization of the operator —dg with periodic boundary conditions. It has the eigenvectors 



ikei _ ( JkBi 



[e"^-«] = {e 



ik9r 



y 



(3.14) 



with entries given by the restriction of the continuum eigenfunctions e'*^^ at the primary grid angles. Here 
k is integer, satisfying |fc| < {n— l)/2, and the eigenvalues are w^, where 



UJk = |fc| 



smc 



fcft.fi 



(3.15) 



and sinc(a;) = sin(a;)/x. Note that uj^ « fc^ only for |fc| <g; n. 

To determine the spectral decomposition of the discrete DtN map A^ we proceed as in the continuum 
and write the potential \Jj as a Fourier sum 



U, =«b(0)1b+ Y. V,{k)[e''''], 

|fe|<^,fe^O 



(3.16) 



where we recall that le G K" is a vector of all ones. We obtain the finite difference equation for the 
coefficients Vj{k), 



1 (V,+,{k)-Vj{k) V,{k)-V,^,{k) 



ai 



"7-1 



-i.lV,ik)=0, 



where j = 2,3, . . . ,£. It is the discretization of ( |3.2[ ) that takes the form 

d f dv{z,ky 



dz \ dz 



-k''v{z,k) = 0, 



(3.17) 



(3.18) 



in the coordinates (3.9), where we let in an abuse of notation v{r, k) -^ v{z, k). The boundary condition at 



r = is mapped to 



lim v{z, k) = 0, 



(3.19) 



and it is implemented in the discretization as Ve+i{k) = 0. At the boundary r = 1, where z — 0, we specify 
Vi{k) as some approximation of vslk). 

The discrete DtN map A^ is diagonahzed in the basis {[e*'^^]}|j.|<ri-i , and we denote its eigenvalues by 



2.2 



F{ujf.). Its definition depends on the type of grid that we use, indexed by fni/2, as explained in section 
In the case mi/2 — 0, the first radius next to the boundary is r2, and we define the boundary flux at r i = 1 
as {Vi{k) — V2{k)) / ai. When mi 12 — 1, the first radius next to the boundary is r2, so to compute the fiux 



at Ti we introduce a ghost layer at j'g > 1 and use equation (3.171 for j = 1 to define the boundary flux as 

VQ{k)-Vi{k) 



a„ 



aiLoiVi{k) 



Vi{k)-V2{k) 



Therefore, the eigenvalues of the discrete DtN map are 



p, 2^ ^ 2 , Vl{k)~V2{k) 



(3.20) 



3.1.2 The NtD map 

The NtD map Aj^ has eigenfunctions e''^^ for /c ^ and eigenvalues f^{k^) — l//(fc^). Equivalently, in terms 
of the solution v{z, k) of equation (3.18) with boundary conditions ( |3.19| and 

dv{0, k) ] '■^'^ 



dz 



we have 



27r 



fik') 



M9)e-'^"d9 = ^^ik), 



v{0,k) 



ifi^ik) 



(3.21) 



(3.22) 



In the discrete case, let us use the grids with mi/2 — 1- We obtain that the potential Vj{k) satisfies (3.17) 
for j = 1,2,...,^, with boundary conditions 



Vi{k)^Vo{k) 



"0 



= <i>^(fc), y,+i=o. 



(3.23) 



Here ^sik) is some approximation of (pig(k). The eigenvalues of Aj^ are 



fH^D 



Vijk) 
$,(fc)- 



(3.24) 



3.2 Rational approximations, optimal grids and reconstruction mappings 

Let us define by analogy to (3.22) and (3.24[) the functions 



ffn^-^M 



fW = 



fe 






(3.25) 



10 



where v solves equation (3.181 with fc^ replaced by A and Vj solves equation (3.171 with cj^ replaced by A 



The spectral parameter A may be complex, satisfying A G C \ (— oo,0]. For simplicity, we suppress in the 
notation the dependence of v and Vj on A. We consider in detail the discretizations on grids indexed by 
"I1/2 = Ij but the results can be extended to the other type of grids, indexed by mi/2 — 0. 



Lemma 1. The function /^(A) is of form 



fw = 



-,^ A — i 



(3.26) 



where fi{t) is the positive spectral measure on (— cxd,0] of the differential operator d-^dz, with homogeneous 



Neumann condition at z = and limit condition (3.19). The function F^ {X) has a similar form 



Ft(A)= / 

J — c 



\-t 



(3.27) 



where /i (i) is the spectral measure of the difference operator in (3.11) with boundary conditions (3.23) 



Proof: The result (3.26) is shown in 03] and it says that /^(A) is essentially a Stieltjes function. To 
derive the representation (3.27), we write our difference equations in matrix form for V ~ [Vi, . . . , V^)"^, 



(A - AI) V 



^h(A ). 
Si 



-ei. 



(3.28) 



Here I is the i x i identity matrix, ei = (1, . . . , 0)"^ G M^ and A is the tridiagonal matrix with entries 



A, 



Oii \ Q.i 



a,-l / ^*.J 



'.J—Si , + ^J—52 n 



=^^^(5,_i , + .^(5^+1 ,- if 1 < i < ^, 1 < 7 < £, 

if ?: = 1, !<]<(.. 



(3.29) 



The Kronecker delta symbol 5ij is one when i = j and zero otherwise. Note that A is a Jacobi matrix when 
it is defined on the vector space M with weighted inner product 



(a, b) = 2^ ^j'^jbjj ^ = {o-ii 



That is to say, 



A = diag(aj' ,...,a/ jAdiagfS^ ,.. 



ibi,...,bi)^ 



-1/2 



(3.30) 



(3.31) 



is a symmetric, tridiagonal matrix, with negative entries on its diagonal and positive entries on its up- 
per/lower diagonal. It follows from [3^ that A has simple, negative eigenvalues —6"^ and eigenvectors 



Yj = (Yij, . . . , Yij) that are orthogonal with respect to the inner product (3.30). We order the eigenval- 
ues as 

Si<62<...<Si, (3.32) 
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and normalize the eigenvectors 



p=l 



(3.33) 



Then, we obtain from (3.25) and (3.28), after expanding V in the basis of the eigenvectors, that 






This is precisely (3.271, for the discrete spectral measure 



(3.34) 



M^(t) 



E 



^,H{-t-d^), i 



where H is the Heaviside step function. D. 



y2 



(3.35) 



Note that any function of the form (3.34) defines the eigenvalues F^(a;^) of the NtD map Aj^ of 



finite volumes scheme with £ + 1 primary radii and uniform discretization in angle. This follows from the 
decomposition in section [3. 1| and the results in |44j . Note also that there is an explicit, continued fraction 
representation of F^{X), in terms of the network conductances, i.e., the parameters aj and Sj, 



Ft(A) = 



(3.36) 



SiA + 



ai 



aeX 



at 



This representation is known in the theory of rational function approximations [591 I44j and its derivation is 
given in appendix [B) 

Since both P{X) and i^^(A) are Stieltjes functions, we can design finite volume schemes (i.e., layered net- 
works) with accurate, rational approximations F\X) of /^(A). There are various approximants i^^(A), with 
different rates of convergence to /^(A), as £ — >■ oo. We discuss two choices below, in sections 3.2.2 



and 



3.2.3 



but refer the reader to [301 UHl 132] for details on various Fade approximants and the resulting discretization 
schemes. No matter which approximant we choose, we can compute the network conductances, i.e., the 
parameters aj and cij for j = !,...,£, from 2£ measurements of /' (A). The type of measurements dictates 
the type of approximant, and only some of them are directly accessible in the EIT problem. For example, 
the spectral measure /i(A) cannot be determined in a stable manner in EIT. However, we can measure the 
eigenvalues /^(fc^) for integer k, and thus we can design a rational, multi-point Fade approximant. 



Remark 1. We describe in detail in appendixUAhow to determine the parameters {ofj, aj}j=i,.../ from 2£ 
point measurements of p{X), such as f^{k'^), for k — 1,.. ., ^^^ — 2£. The are two steps. The first is to 
write f^(A) as the ratio of two polynomials of X, and determine the 2£ coefficients of these polynomials from 



the measurements F^{ijj^) of f^{k'^), for 1 < k < ^^^- See section 3.2.2 for examples of such measurements. 



12 



The exponential instability of EIT conies into play in this step, because it involves the inversion of a Van- 
dermonde matrix. It is known I33}j that such matrices have condition numbers that grow exponentially with 
the dimension L The second step is to determine the parameters {a.j,aj}j=i,,,,^( from the coefficients of the 
polynomials. This can he done in a stable manner with the Euclidean division algorithm JJj 



The approximation problem can also be formulated in terms of the DtN map, with F{\) = 1/F^{X). 



Moreover, the representation ( 3.36 1 generalizes to both types of grids, by replacing Si A with aiTOi/2A. Recall 
equation (3.201 and note the parameter Si does not play any role when mi/2 — 0. 



3.2.1 Optimal grids and reconstruction mappings 

Once we have determined the network conductances, that is the coefficients 



dr 



r,+i ra{r) '' 



"^^ (j{r) 



dr. 



J ^ !,...,£, 



(3.37) 



'•3+1 



we could determine the optimal placement of the radii rj and fj, if we knew the conductivity a{r). But a{r) 
is the unknown in the inverse problem. The key idea behind the resistor network approach to inversion is 
that the grid depends only weakly on a, and we can compute it approximately for the reference conductivity 

a(°) = 1. 



Let us denote by /^'■"■'(A) the analog of (3.25) for conductivity ct*^°\ and let F^^°^{X) be its rational 
approximant defined by (3.36), with coefficients a° and a" given by 



^J + 1 



i(°) - 



(°) r 



' dr , r) ' 
— = log -^-—. 

r Wo) 



3 = h- 



(3.38) 



-(°) 



Since r^ = r-j^ = 1, we obtain 



(o) 



exp 



9=1 



„(°) 






exp 



E< 

9=1 



:(°) 



J = l, 



(3.39) 



We call the radii (3.39) optimal. The name refers to the fact that finite volume discretizations on grids with 



such radii give an NtD map that matches the measurements of the continuum map A („, for the reference 
conductivity a''°\ 

Remark 2. It is essential that the parameters {aj,aj} and {a ° ,a ° } are computed from the same type of 
measurements. For example, if we measure f^{k'^), we compute {aj, Sj} so that 

FHcol) = fHk% 

and {a " ,a° } so that 



where k — 1, . . . ,{n — l)/2. This is because the distribution of the radii {3.39) in the interval [0, 1] depends 
on what measurements we make, as illustrated with examples in sections\3.2.S\ and\3.2.3{ 
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Now let us denote by I>„ the set in M 2 of measurements of /^(A), and introduce the reconstruction 

71-1 

mapping Q„ defined on I?„, with values in M^^ . It takes the measurements of /^(A) and returns the 
(n — l)/2 positive numbers 



7" + l-mi/2 ~ 




j = 2 - mi/2, 


^J+mi/2 - 




j = l,2,...,£ 



(3.40) 



where we recall the relation (2.14) between i and n. We call Q„ a reconstruction mapping because if we 



r* ^ ' o VI H t^ ' 



take aj and CTj as point values of a conductivity at nodes r, and f;- , and interpolate them on the optimal 
grid, we expect to get a conductivity that is close to the interpolation of the true cr(r). This is assuming that 
the grid does not depend strongly on (7{r). The proof that the resulting sequence of conductivity functions 
indexed by i converges to the true a{r) as ^ — > 00 is carried out in [T3], given the spectral measure of /^(A). 
We review it in section |3.3[ and discuss the measurements in section |3.2.3| The convergence proof for other 
measurements remains an open question, but the numerical results indicate that the result should hold. 
Moreover, the ideas extend to the two dimensional case, as explained in detail in sections |4] and [5] 

3.2.2 Examples of rational interpolation grids 

Let us begin with an example that arises in the discretization of the problem with lumped current measure- 
ments 



1 fl^q+i 

J, = — / A„UB{e)de, 



for he — ^^, and vector Ug = (wb(0i), . . . , UB{dn)) of boundary potentials. If we take harmonic boundary 
excitations ub{0) = e'*^^, the eigenfunction of A^ for eigenvalue fik"^), we obtain 



J, = Y L '^' Ke'^'de = f{k^) 



khfi 
sine , 

2 



ikOa _ f\'^ ), , ike„ 



e'""' = '-^jj^coke^'''- , q = l,...,n. (3.41) 

\k\ 



These measurements, for all integers k satisfying |fc| < ^^y^, define a discrete DtN map M„(Act). It is a 
symmetric matrix with eigenvectors [e'*^^] = (e*'^^^, . . . , e**^^") , and eigenvalues -^L', ' ujk- 

The approximation problem is to find the finite volume discretization with DtN map A-y = M„(Acr). 
Since both A^ and M„ have the same eigenvectors, this is equivalent to the rational approximation problem 



of finding the network conductances (3.8 1 (i.e., aj and aj), so that 



.2^_/(fc'),, ._. n-1 



^K) = ^^fc, k^l,...,^. (3.42) 

The eigenvalues depend only on |fc|, and the case fc = gives no information, because it corresponds to 
constant boundary potentials that lie in the null space of the DtN map. This is why we take in (|3.42 1 only 



the positive values of fc, and obtain the same number (n — l)/2 of measurements as unknowns: {aj}j=i,...,f 
and {Sj}j=2-mi/2, ...,£• 

When we compute the optimal grid, we take the reference cr^°-' = 1, in which case /'"'(fc^) = |fc|. Thus, 
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m^5, m ^1, n^25 



o X o X o xo* 
o X o X o xc* 



m=8, m^i^=0, n=35 

X o X o X o X oxoxoxoxc 
X o X o X o X oxoxoxoxc 



0.4 0.6 

re [0,1] 



1 



0.2 



0.4 0.6 

re [0,1] 



0.8 



Figure 3: Examples of optimal grids with n equidistant boundary points and primary and dual radii shown 
with X and o. On the left we have rt = 25 and a grid indexed by ?7ii /2 = 1, with £ = m + 1 = 6. On the right 
we have n — 35 and a grid indexed by ?7ii/2 = 0, with i = m + 1 = 8. The grid shown in red is computed 
with formulas (3.44). The grid shown in blue is obtained from the rational approximation (3.50). 



the optimal grid computation reduces to that of rational interpolation of /(A), 

F(°)(c.D=c., = /(°)(c.D, ^ = 1,...,^. 



(3.43) 



This is solved explicitly in [TU]. For example, when mi/2 = 1, the coefficients a° and a° are given by 



„(°) 



hg cot 



-i2£~2j + l) 



a J = hg cot 



-i2£~2j + 2) 



and the radii follow from (3.39). They satisfy the interlacing relations 



1 = f<°) = r[°'> > f<°) > 4°) > . . . > f<;\ > r^;\ > 0, 



J ^1,2. ..J, (3.44) 



(3.45) 



as can be shown easily using the monotonicity of the cotangent and exponential functions. We show an 
illustration of the resulting grids in red, in Figure [3] Note the refinement toward the boundary r = 1 and 
the coarsening toward the center r = of the disk. Note also that the dual points shown with o are almost 
half way between the primary points shown with x . The last primary radii r^Vi are small, but the points 
do not reach the center of the domain at r = 0. 

In sections [4] and [5] we work with slightly different measurements of the DtN map A-y = M„(Acr), with 
entries defined by 






(3.46) 



using the non-negative measurement (electrode) functions Xq{(^) that are compactly supported in {9q,9q+i)^ 
and arc normalized by 



For example, we can take 



xq{e)de^i. 



I 0, otherwise. 



15 



and obtain after a calculation given in appendix [O that the entries of A^ are given by 

' khfi 



{A,)^,=lj:^''^'''''^fik') 



kez 



P,9,= 1, 



We also show in appendix [C] that 



A. [e'"'] = ^FU) [e'^1 , \k\ < ^.. 



with eigenvectors [e*'^^] defined in (3.14) and scaled eigenvalues 



Fiu^l) = f{k^) 



sine 



khf. 



n 2 



= fU) 



sine 



kha 



Here we recalled (3.42 1 and (3.15) 



There is no explicit formula for the optimal grid satisfying 

i?(°)(..^)^F('')(.^)sinc('^') =... 



khf, 



(3.47) 



(3.48) 



(3.49) 



(3.50) 



but we can compute it as explained in Remark [T] and appendix [D] We show in Figure [3] two examples of the 
grids, and note that they are very close to those obtained from the rational interpolation (|3.43 ) . This is not 



surprising because the sine factor in (3.50) is not significantly different from 1 over the range |fc| < ^^^, 



2 , «in[f(l-^)] ^ 



smc 



khf. 



< 1. 



Thus, many eigenvalues F^°'{uyf,) are approximately equal to Wfc, and this is why the grids are similar. 



3.2.3 Truncated measure and optimal grids 

Another example of rational approximation arises in a modified problem, where the positive spectral measure 
jjL in Lemma [l] is discrete 

oo 

m^~Y.^,H{-t-5^). (3.51) 



i=i 



This does not hold for equation (^3.2| or equivalently (3.18), where the origin of the disc r = is mapped to 



cxD in the logarithmic coordinates z{r)^ and the measure /x(i) is continuous. To obtain a measure like (3.51 ), 
we change the problem here and in the next section to 



a{r) dr 
with e e (0, 1) and boundary conditions 



r(j(r) 



dv{r) 
dr 



A«(r)=0, re(e,l). 



dv{o) 
dr 



Vb, w(e) = 0. 



(3.52) 



(3.53) 
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The Dirichlet boundary condition at r = e may be realized if we have a perfectly conducting medium in the 
disk concentric with fl and of radius e. Otherwise, v{e) — gives an approximation of our problem, for small 
but finite e. 

Coordinate change and scaling 

It is convenient here and in the next section to introduce the scaled logarithmic coordinate 

CM = ^ = I /' f ' ^ = -l«g(^) = ^^''^(^)' (3-54) 



and write (3.9) in the scaled form 



'^'^ - f' ^' -z'(C), M= /^(,(t))dt^z'(C). (3.55) 



Z J, a{r{t)) '^'■- Z 

The conductivity function in the transformed coordinates is 

<^'{C) = <y{r{C)), KC) = e-^<, (3.56) 

and the potential 

V (z) = - 

satisfies the scaled equations 

'^"(^^ - 1, .(L') = 0, (3.58) 



v'iz') = ^^M^ (3.57) 



dz 
where we let A' = X/Z'^ and 



"^ dt 
Remark 3. We assume in the remainder of this section and in section \3.S\ that we work with the scaled 



L' = z'{l)= f -^ (3.59) 

Jo <^ (t) 



equations {3.58) and drop the primes for simplicity of notation. 



The inverse spectral problem 

The differential operator ^ -^ acting on the vector space of functions with homogeneous Neumann conditions 
at z = and Dirichlet conditions at z = L is symmetric with respect to the weighted inner product 

(a, 6)= [ a{z)h{z)dz^ f a{z{C))b{z{C))cr{OdC, L = z(l). (3.60) 

Jo Jo 
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It has negative eigenvalues {— (5^}j=i.2,..., the points of increase of the measure (3.511, and eigenfunctions 



yj{z). They are orthogonal with respect to the inner product (3.601, and we normalize them by 



W^l 



iyj-.V]) 



y^{z)dz^ 1. 



The weights S,j in (3.51) are defined by 



0=2/|(0). 



(3.61) 



(3.62) 



For the discrete problem we assume in the remainder of the section that mi/2 = 1, and work with the 
NtD map, that is with F^{X) represented in Lemma ll] in terms of the discrete measure fJ,^{t). Comparing 



( 3.51 ) and ( 3.35 ) , we note that we ask that /i^ (i) be the truncated version of fi{t) , given the first £ weights ^j 
and eigenvalues ~S'^, ior j — 1, . . . , I. We arrived at the classic inverse spectral problem [MJ [121 IMl [M1I55) . 
that seeks an approximation of the conductivity a from the truncated measure. We can solve it using the 
theory of resistor networks, via an inverse eigenvalue problem |20j for the Jacobi like matrix A defined in 



(3.29). The key ingredient in the connection between the continuous and discrete eigenvalue problems is the 



optimal grid, as was first noted in [I^] and proved in ^14]. We review this result in section 3.3 



The truncated measure optimal grid 



The optimal grid is obtained by solving the discrete inverse problem with spectral data for the reference 
conductivity a'^°\(), 

2?i°) = [^f ^2,sf=n(^j-\y J - 1, . . . , fj . (3.63) 



,» ;;;(°)l 



The parameters {a^ , (ij}j=i^,,,^e can be determined from I?„ with the Lanczos algorithm [65l[20] reviewed 
briefly in appendix IE] The grid points are given by 



q=l 9=1 



J = l 



f 

, . . . , -o, 



(3.64) 



where q" = Q° — 0. This is in the logarithmic coordinates that are related to the optimal radii as in 



(3.56). The grid is calculated explicitly in [14, Appendix A]. We summarize its properties in the next lemma, 
for large £. 

Lemma 2. The steps {a°,a°}j=i^,,,,i of the truncated measure optimal grid satisfy the monotone relation 



<o) ^ (o) ^ ^(o) ^ (o) 

Moreover, for large i, the primary grid steps are 



a 



(o) 



V2+o{e-^) 



<o) ^ (o) 



(3.65) 



(3.66) 
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Figure 4: Example of truncated measure optimal grid with £ = 6. This is in the logarithmic scaled 
coordinates C G [0, 1]. The primary points are denoted with x and the dual ones with o. 



xOxO xO xO xO xO 

xOxOxO xO xO xO 

)0 kO xO X O X O X O 







0.2 



0.4 



0.6 



0.8 



Figure 5: The radial grid obtained with the coordinate change r — e~^^ . The scale Z = — loge affects 
the distribution of the radii. The choice e = 0.1 is in blue, e — 0.05 is in red and e = 0.01 is in black. The 
primary radii are indicated with x and the dual ones with o. 



and the dual grid steps are 



i-j)-'+r'] 



j(o) ^ 2 + o [je 

vrv/^2 _ ^j _ 1/2)2 



1 < J < ^• 



(3.67) 



We show in Figure |4] an example for the case i — 6. To compare it with the grid in Figure [3j we plot in 
Figurep^the radii given by the coordinate transformation (3.56), for three different parameters e. Note that 
the primary and dual points are interlaced, but the dual points are not half way between the primary points, 
as was the case in Figure [3] Moreover, the grid is not refined near the boundary at r = 1. In fact, there is 
accumulation of the grid points near the center of the disk, where we truncate the domain. The smaller the 
truncation radius e, the larger the scale Z = — log e, and the more accumulation near the center. 

Intuitively, we can say that the grids in Figure [3] are much superior to the ones computed from the 
truncated measure, for both the forward and inverse EIT problem. Indeed, for the forward problem, the rate 
of convergence of F^{X) to p{X) on the truncated measure grids is algebraic [T3] 



|/t(A)-Ft(A)| = 



E 







A + 52 




The rational interpolation grids described in section 3.2.2 give exponential convergence of F^ {X) to /^(A) 



[51j . For the inverse problem, we expect that the resolution of reconstructions of a decreases rapidly away 
from the boundary where we make the measurements, so it makes sense to invert on grids like those in Figure 
[31 that are refined near r = \. 

The examples in Figures [3] and [S] show the strong dependence of the grids on the measurement setup. 
Although the grids in Figure [5] are not good for the EIT problem, they are optimal for the inverse spectral 
problem. The optimality is in the sense that the grids give an exact match of the spectral measurements 
(3.63) of the NtD map for conductivity a'^°\ Furthermore, they give a very good match of the spectral 
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measurements (3.681 for the unknown a, and the reconstructed conductivity on them converges to the true 



a, as we show next. 

3.3 Continuum limit of the discrete inverse spectral problem on optimal grids 

Let Qn : Vn — > M^^ be tlie reconstruction mapping that takes the data 

Vn^{i„5,, i = !,...,£} (3.68) 

to the 211 — ^^yi positive values {cj, (Tj}j=i,...,£ given by 



a, ^ a 



(o) 



(7, = ^, dj+i = ^—, j = l,2,...,i. (3.69) 



The computation of {oij,o.j}j=i,...,i requires solving the discrete inverse spectral problem with data T>n, 
using for example the Lanczos algorithm reviewed in appendix [El We define the reconstruction cr^{C) of the 



conductivity as the piecewise constant interpolation of the point values (3.691 on the optimal grid (3.64 1. 
We have 

ifce[ci°\g:-\), j-i,...,^, 

^'(C) = <( a„ ifCe[g°\CJ°^), J=2,...,£+l, (3.70) 

^e+1, ifCe[cS,i] 

and we discuss here its convergence to the true conductivity function ^{Q, as -^ — > oo. 

To state the convergence result, we need some assumptions on the decay with j of the perturbations of 
the spectral data 

A6,^S,^6^''\ A^,=^,-^^\ (3.71) 

The asymptotic behavior of Sj and ^j is well known, under various smoothness requirements on a{z) [5511601 
W] . For example, if o-(C) G H^[0, 1], we have 

A<5, = 5, - Sl"^ = j^0^ + O ir') and A^, = ?, - C]"^ - O (,-) , (3.72) 



where g(C) is the Schrodinger potential 



<z(C)=a(C)-^. (3.73) 



We have the following convergence result proved in [H] . 

Theorem 1. Suppose that cr(C) is a positive and bounded scalar conductivity function, with spectral data 
satisfying the asymptotic behavior 



.-loW)- "''^"(j- 



A6j = O ^- — —r I , A^j = O ( — ) , for some s > 1, as j ^ oo. (3.74) 
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Then cr^(C) converges to cr(C) as £ ^ oo, pointwise and in L^[Q, 1]. 



Before we describe the outline of the proof in [14], let us note that it appears from (3.72) and (3.74) that 



the convergence result applies only to the class of conductivities with zero mean potential. However, if 



q{C)dC + 0, 



(3.75) 



we can modify the point values (3.69) of the reconstruction u (C) by replacing or° and a° with a 



(o) 



;(o) 



„(«) 



and a,' , for j = 1, . . . ,^. These are computed by solving the discrete inverse spectral problem with data 



^-{^.^ 






■(«) 



j = i, 



, , £ > , for conductivity function 



-^^nc)4( 



e^C + g-y«C 



This conductivity satisfies the initial value problem 



rf^y^^io 



dC^ 



and we assume that 



Wf7(9)(C) for 0<C<1, 



dcr(9)(0) 
dC, 



and o-(«)(0) = l 



so that (3.76) stays positive for Q e [0, 1]. 



9> 



r(9) 



t(9) 



(5 and ^j — o satisfy the assumptions (|3.74|), so Theorem 



(3.76) 



(3.77) 



(3.78) 



As seen from (3.72), the perturbations 8j 
[ij applies to reconstructions on the grid given by ct^^'. We show below in Corollary [ij that this grid is 
asymptotically the same as the optimal grid, calculated for cr(°\ Thus, the convergence result in Theorem [ij 



applies after all, without changing the definition of the reconstruction (3.70) 



3.3.1 The case of constant Schrodinger potential 

The equation (3.58) for a -^ cr® can be transformed to Schrodinger form with constant potential q 

d^wiC) 



dC 



dw{0) 



0, CG(0,1), 



(3.79) 



= -1, 



u;(l)=0. 



by letting w(^) = w(C)-\/cr(9)(^). Thus, the eigenfunctions y^ (C) of the differential operator associated with 



cr^'^-'(C) are related to y° {(,), the eigenfunctions for a'^°^ = 1, by 






2/1°^ (C) 



®(C) 



(3.80) 



They satisfy the orthonormality condition 



' yf (C)y® (C)<T® (C)rfC = /' yf\0y^^\0dC = <5.p, 

Ja 



(3.81) 



21 



and since (t('')(0) = 1, 



e 



(5) 



m. 



y^'iO) = y]"^(0) 



,(o)/ 



&\ J = 1,2,.. 



(3.82) 



The eigenvalues are shifted by q, 



C®^' = _f^(o)'' 



J = l,2,.. 



(3.83) 



Let {aij^ ,a '^ }j=i^,,,^i be the parameters obtained by solving the discrete inverse spectral problem with 



data Vn' . The reconstruction mapping Q„ : Vn' 



P^ gives the sequence of 2£ = ^^ pointwise values 



M _ "j 



(5) 



(o) 



,^. 



(3.84) 



We have the following result stated and proved in [T^ . See the review of the proof in appendix [F| 

Lemma 3. The point values ay' satisfy the finite difference discretization of initial value problem (3.77), 
on the optimal grid, 



a 



1 

Jo) 



/9) 

i+1 



J-1 



a 



(o) 



^(°^ 

^i-1 



-qjaf = 0, J =2, 3 



, ... J -I., 



<y^ — \ a. 



^i") 



„(°) 



qy^?' 



0, af) = 1. 



(3.85) 



Moreover, af_l, = ^'^f'^f+i, for j = l,..., I 

The convergence of the reconstruction a^''>'^{C,) follows from this lemma and a standard finite-difference 
error analysis [36] on the optimal grid satisfying Lemmal2] The reconstruction is defined as in (3.70), by the 
piecewise constant interpolation of the point values ( 3.84[ ) on the optimal grid. 

Theorem 2. j4s ^ ^- oo we have 



max 



af -a®(CJ°^) 



and max 



;(?) 



;?]^i-^^'^(0:i) 



{q)(7(°) 



0, 



(3.86) 



and the reconstruction a^'' >'{(,) converges to cr^'''{C,) *^ L°°[Q, 1]. 

As a corollary to this theorem, we can now obtain that the grid induced by cr®(C), with primary nodes 
C and dual nodes Cf, is asymptotically close to the optimal grid. The proof is in appendix F 

Corollary 1. The grid induced by cr^^'{C) is defined by equations 



^®(C) ^. " ' 7o ^^^ ^ h " 



J = !,...,£, Cf)=cj^^=0, (3.87) 



p=i 



p=i 



and satisfies 



max 



cf cj°' 



0, 



max 

i<j<e+i 



7(q) _ 7{o) 



0, as ^ —> oo. 



(3.88) 
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3.3.2 Outline of the proof of Theorem [T] 



The proof given in detail in |14| has two main steps. The first step is to establish the compactness of the 
set of reconstructed conductivities. The second step uses the established compactness and the uniqueness of 
solution of the continuum inverse spectral problem to get the convergence result. 

Step 1: Compactness 



We show here that the sequence {c (C)}f>i of reconstructions (3.701 has bounded variation 



Lemma 4. The sequence {cTj, (Tj+i}j=i..../ [3.69) returned by the reconstruction mapping Q„ satisfies 



^ \\ogdj+i - logo-jl + ^ |log%+i - logcTj+il < C, 



(3.89) 



where C is independent of i. Therefore the sequence of reconstructions {c (C)}^>i has uniformly hounded 
variation. 



Our original formulation is not convenient for proving (3.891, because when written in Schrodinger form, 



it involves the second derivative of the conductivity as seen from (|3.73[). Thus, we rewrite the problem in 



first order system form, which involves only the first derivative of o'(C), which is all we need to show (3.891 
At the discrete level, the linear system of (. equations 



AV- AV 



6]^ 

Si 



\T ■ 



for the potential V — {Vi, . . . , Vi) is transformed to the system of 2£ equations 

BH5W \/AH3W= ^^ 



'Xai 



for the vector W = [Wi,W2 



(Wi,W2, . . . ,We,Wi+iy 



W, 



/^V,, W,+i 



with components 



^J + l ( ^J + l - ^3 



/Xa-j 



v(°) 



, J = l, 



Here H = diag I S^ , cti , . . . , a/, a\° j and B is the tridiagonal, skew-symmetric matrix 



B 






/3i 








-/3i 





/32 








-/32 








'f32e-i J 



(3.90) 



(3.91) 



(3.92) 



(3.93) 
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with entries 



J2p 



v/«P«P+i ^4°)a^°|i V '^P 



CTp+l _ (o) (Tp+1 
^P V CTp+l ■ 



ft 



'2p-l 



Cp+l _ (o) / gp+1 



Note that we have 

2£-l 



E 



log 



4°^ 



^^^ y^^v '"^ 



1 1 



(3.94) 
(3.95) 



(3.96) 



p=i 



p=i 



and we can prove (3.89) by using a method of small perturbations. Recall definitions (3.71 1 and let 



AS^^^rAS,, Ae;: = rA6, j = 1, 



(3.97) 



where r € [0, 1] is an arbitrary continuation parameter. Let also /SJ be the entries of the tridiagonal, 
skew-symmetric matrix B'' determined by the spectral data S^ = 5° + A(5^ and ^'' = C° + A^'', for 
i — !,...,£. We explain in appendix [Gl how to obtain explicit formulae for the perturbations d\ogf3^ in 
terms of the eigenvalues and eigenvectors of matrix B'' and perturbations d6^ = A6jdr and d^^ = A^jdr. 
These perturbations satisfy the uniform bound 



2i-l 

Y^ |dlog/3;| <Ci\dr 



(3.98) 



with constant Ci independent of £ and r. Then, 

log 



/?. 



^r^ ■ - 



dlog/3^ 



(3.99) 



2«-l 



satisfies the uniform bound 2, 
Step 2: Convergence 



log 






< Ci and (3.891 follows from (3.961. 



Recall section 3.2 where we state that the eigenvectors Y^ of A are orthonormal with respect to the weighted 
inner product ( 3.30 ) . Then, the matrix Y with columns diag I a j^ , . . . , a| ) Y^ is orthogonal and we have 



(YY^)n-"iE^. = l- 



(3.100) 



j=i 
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Similarly 



4°^Ee 



(o) _ 



2l'i 






= 1, 



J=l 



where we used (3.63), and since A^j are summable by assumption (3.74), 



(3.101) 



CTi 






S(°) 



i°^E^^. 



1 + 0(3^°^) = 1 + 



i=i 



(3.102) 



But cr^(O) = (Ti, and since cr^(C) has bounded variation by Lemma kl we conclude that cr^(C) is uniformly 
bounded in C G [0, 1]. 

Now, to show that (t^(C) — )■ cr(C) in i^[0, 1], suppose for contradiction that it does not. Then, there exists 
£ > and a subsequence a^'' such that 

Ik*"" -ct||li[o,i] ^ ^■ 

But since this subsequence is bounded and has bounded variation, we conclude from Helly's selection principle 
and the compactness of the embedding of the space of functions of bounded variation in L^ [0, 1] [57] that it 
has a convergent subsequence pointwise and in L^[0, 1]. Call again this subsequence a^'' and denote its limit 
by a* 7^ a. Since the limit is in i^[0, 1], we have by definitions (3.55) and Remark^ 



z(C;a^ 



dt 



<^''{t) 



<C;<^) 



dt 
^W)' 



z{Q-a'^) 



(7^''it)dt^z{C 



a*) ^ I a{t)dt. (3.103) 

"'0 



Furthermore, the continuity of p with respect to the conductivity gives /^(A;(T^*) — >■ /^(A;cr*). However, 
Lemmallland (3.51) show that /^(A;cr^) -^ p{X;a) by construction, and since the inverse spectral problem 
has a unique solution [35l |49l [21] |60] , we must have a* — a. We have reached a contradiction, so (J^{() — >■ cr(C) 
in i^[0, 1]. The pointwise convergence can be proved analogously. 



Remark 4. All the elements of the proof presented here, except for establishing the bound (3.98), apply 



to any measurement setup. The challenge in proving convergence of inversion on optimal grids for general 
measurements lies entirely in obtaining sharp stability estimates of the reconstructed sequence with respect to 
perturbations in the data. The inverse spectral problem is stable, and this is why we could establish the bound 



[3.98). The EIT problem is exponentially unstable, and it remains an open problem to show the compactness 



of the function space of reconstruction sequences a from measurements such as {3.49) 



4 Two dimensional media and full boundary measurements 

We now consider the two dimensional EIT problem, where a — a{r,0) and we cannot use separation of 
variables as in section [3| More explicitly, we cannot reduce the inverse problem for resistor networks to 
one of rational approximation of the eigenvalues of the DtN map. We start by reviewing in section |4.1| 
the conditions of unique recovery of a network (r,7) from its DtN map A-y, defined by measurements of 
the continuum A^-. The approximation of the conductivity a from the network conductance function 7 is 
described in section [4?2l 
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4.1 The inverse problem for planar resistor networks 

The unique recoverability from A^ of a network (r,7) with known circular planar graph F is established in 
PSI [^ [231 US]- A graph F = (T', £) is called circular and planar if it can be embedded in the plane with 
no edges crossing and with the boundary nodes lying on a circle. We call by association the networks with 
such graphs circular planar. The recoverability result states that if the data is consistent and the graph F is 
critical then the DtN map A^ determines uniquely the conductance function 7. By consistent data we mean 
that the measured matrix A-y belongs to the set of DtN maps of circular planar resistor networks. 

A graph is critical if and only if it is well- connected and the removal of any edge breaks the well- 
connectedness. A graph is well-connected if all its circular pairs {P, Q) are connected. Let P and Q be two 
sets of boundary nodes with the same cardinality \P\ — \Q\. We say that (P, Q) is a circular pair when the 
nodes in P and Q lie on disjoint segments of the boundary B. The pair is connected if there are \P\ disjoint 
paths joining the nodes of P to the nodes of Q. 

A symmetric n x n real matrix A^ is the DtN map of a circular planar resistor network with n boundary 
nodes if its rows sum to zero A-yl = (conservation of currents) and all its circular minors (A^)pq have 
non-positive determinant. A circular minor (A^)p.q is a square submatrix of A^ defined for a circular 
pair (P, Q), with row and column indices corresponding to the nodes in P and Q, ordered according to a 
predetermined orientation of the circle B. Since subsets of P and Q with the same cardinality also form 
circular pairs, the determinantal inequalities are equivalent to requiring that all circular minors be totally 
non-positive. A matrix is totally non-positive if all its minors have non-positive determinant. 



Examples of critical networks were given in section 2.2 with graphs F determined by tensor product 
grids. Criticality of such networks is proved in |22j for an odd number n of boundary points. As explained 
in section 2.2 (see in particular equation (2.14)), criticality holds when the number of edges in £ is equal to 
the number n{n — l)/2 of independent entries of the DtN map A^. 

The discussion in this section is limited to the tensor product topology, which is natural for the full 
boundary measurement setup. Two other topologies admitting critical networks (pyramidal and two-sided), 
are discussed in more detail in sections |5.2.1| and |5.2.2[ They are better suited for partial boundary mea- 
surements setups [m [17] . 



Remark 5. It is impossible to recover both the topology and the conductances from the DtN map of a 
network. An example of this indetermination is the so-called Y — A transformation given in figure [M A 
critical network can be transformed into another by a sequence of Y ~ IS. transformations without affecting 
the DtN map i23j . 

4.1.1 From the continuum to the discrete DtN map 

Ingerman and Morrow [42' showed that pointwise values of the kernel of K^ at any n distinct nodes on 
B define a matrix that is consistent with the DtN map of a circular planar resistor network, as defined 
above. We consider a generalization of these measurements, taken with electrode functions Xg(^)j ^^ given 



in equation (3.46). It is shown in [T3I that the measurement operator M„ in (3.46) gives a matrix M„(Act) 



that belongs to the set of DtN maps of circular planar resistor networks. We can equate therefore 

M„(A,) = A^, (4.1) 
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Figure 6: Given some conductances in the Y network, there is a choice of conductances in the A network 
for which the two networks are indistinguishable from electrical measurements at the nodes p, q and r. 

and solve the inverse problem for the network (F, 7) to determine the conductance 7 from the data A-y. 



4.2 Solving the 2D problem with optimal grids 

To approximate (t(x) from the network conductance 7 we modify the reconstruction mapping introduced in 
section [3^ for layered media. The approximation is obtained by interpolating the output of the reconstruction 
mapping on the optimal grid computed for the reference a^°' = 1. This grid is described in sections 



2.2 



and 



3.2.2 But which interpolation should we take? If we could have grids with as many points as we wish, the 
choice of the interpolation would not matter. This was the case in section [373l where we studied the continuum 
limit n — )• 00 for the inverse spectral problem. The EIT problem is exponentially unstable and the whole idea 
of our approach is to have a sparse parametrization of the unknown a. Thus, n is typically small, and the 
approximation of a should go beyond ad-hoc interpolations of the parameters returned by the reconstruction 
mapping. We show in section |4.2.3| how to approximate a with a Gauss-Newton iteration preconditioned 
with the reconstruction mapping. We also explain briefly how one can introduce prior information about a 
in the inversion method. 



4.2.1 The reconstruction mapping 

The idea behind the reconstruction mapping is to interpret the resistor network (F, 7) determined from the 



measured A^ = M„(Ao.) as a finite volumes discretization of the equation ( 1.1 ) on the optimal grid computed 
for (t'-''' = 1. This is what we did in section 3.2 for the layered case, and the approach extends to the two 
dimensional problem. 

The conductivity is related to the conductances j{E), for E £ £, via quadrature rules that approximate 



the current fiuxes (2.5) through the dual edges. We could use for example the quadrature in [151 [T51 [52] . 
where the conductances are 



7a,h = 0-(Pa,h)- 



(4.2) 



LiEa,b) ' 

{a,b) e {(z,j± i) , (i± i,j')} and L denotes the arc length of the primary and dual edges E and S (see 
section [2.1| for the indexing and edge notation). Another example of quadrature is given in [ 



It is 



specialized to tensor product grids in a disk, and it coincides with the quadrature (3.8) in the case of layered 



media. For inversion purposes, the difference introduced by different quadrature rules is small (see |15[ 
Section 2.4]). 



27 



To define the reconstruction mapping Q„, we solve two inverse problems for resistor networks. One 
with the measured data A^ — M„(Ao.), to determine the conductance 7, and one with the computed data 
A^(o) = M„(Aa ), for the reference cr^°^ = 1. The latter gives the reference conductance 7*-°^ which we 



associate with the geometrical factor in (4.2 1 



so that we can write 



(o) ^(Sg^fc) , . 



<j{Pa,,)^<Ja,,^^. (4.4) 

Tai 



Note that (4.4) becomes (3.40) in the layered case, where (3.8) gives aj — hg/j^^i^^ and Sj — hgj^^^^i. 
The factors hg cancel out. 

Let us call I?„ the set in M*^ of e = n{n — l)/2 independent measurements in M„(Acr), obtained by 
removing the redundant entries. Note that there are e edges in the network, as many as the number of the 
data points in I?„, given for example by the entries in the upper triangular part of M„(Ao-), stacked column 



by column in a vector in M*^. By the consistency of the measurements (section [4.1.1 ), I?„ coincides with 



the set of the strictly upper triangular parts of the DtN maps of circular planar resistor networks with n 
boundary nodes. The mapping Q„ : 2?„ — > M^ associates to the measurements in 2?„ the e positive values 



(TaL in (4.4) 



We illustrate in Figure JTJJb) the output of the mapping Q„, linearly interpolated on the optimal grid. 
It gives a good approximation of the conductivity that is improved further in Figure m c) with the Gauss- 
Newton iteration described below. The results in Figure [7] are obtained by solving the inverse problem for 
the networks with a fast layer peeling algorithm |22] . Optimization can also be used for this purpose, at 
some additional computational cost. In any case, because we have relatively few n(n— l)/2 parameters, the 
cost is negligible compared to that of solving the forward problem on a fine grid. 

4.2.2 The optimal grids and sensitivity functions 

The definition of the tensor product optimal grids considered in sections |2.2| and [3] does not extend to 
partial boundary measurement setups or to non-layered reference conductivity functions. We present here an 
alternative approach to determining the location of the points Pa,b at which we approximate the conductivity 



in the output (4.4) of the reconstruction mapping. This approach extends to arbitrary setups, and it is based 
on the sensitivity analysis of the conductance function 7 to changes in the conductivity |16j . 

The sensitivity grid points are defined as the maxima of the sensitivity functions Dg-^a.bi'^)- They are 
the points at which the conductances 7a, & are most sensitive to changes in the conductivity. The sensitivity 
functions Do-7(x) are obtained by differentiating the identity ^^la) — M„(Ao-) with respect to a, 

(i?<,7)(x)-fi?^A^L_j^ ,, T'vec(M„(i?/C,)(x)), x e ^i. (4.5) 



The left hand side is a vector in W^ . Its fc— th entry is the Frechet derivative of conductance 7^ with respect 
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Figure 7: (a) True conductivity phantoms, (b) The output of the reconstruction mapping Q„, linearly 



interpolated on a grid obtained for layered media as in section 3.2.2 (c) One step of Gauss-Newton improves 
the reconstructions. 



to changes in the conductivity a. The entries of the Jacobian D^A^ e M'^^'^ are 

[D^A 



i^'^i/jk 



dlk 



(4.6) 



where vec(A) denotes the operation of stacking in a vector in W the entries in the strict upper triangular 



part of a matrix A £ M"^". The last factor in (4.5 1 is the sensitivity of the measurements to changes of the 
conductivity, given by 



Xtix)DICa{x; X, y)xj{y)dxdy, 



i^h 



BxB 



(M„(i?/C.))„. (x) = 



^X! / Xr{x)DlC„{yi.]x,y)xk{y)dxdy, i=. 



(4.7) 



Here ICa{x, y) is the kernel of the DtN map evaluated at points x and y £ B. Its Jacobian to changes in the 
conductivity is 



DlC^i^; X, y) = a{x)a{y) {Vx(n(x) • V,G(x, x))} • {Vx(n(y) • V,yG(x, y))} , 



(4.8) 
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Figure 8: Sensitivity functions diag (1/7''^')Z?ct7 computed around the conductivity a — 1 for n = 13. The 
images have a hnear scale from dark blue to dark red spanning ± their maximum in absolute value. Light 
green corresponds to zero. We only display 6 sensitivity functions, the other ones can be obtained by integer 
multiple of 27r/13 rotations. The primary grid is displayed in solid lines and the dual grid in dotted lines. 
The maxima of the sensitivity functions are very close to those of the optimal grid (intersection of solid and 
dotted lines). 

where G is the Green's function of the differential operator u — > V-((TV?i) with Dirichlet boundary conditions, 
and n{x) is the outer unit normal at x e S. For more details on the calculation of the sensitivity functions 
see dH Section 4]. 

The definition of the sensitivity grid points is 



Pa^b — argmax {D„^^ b)(^)j evaluated a,t a — ct^°' = 1. 

X6z0 



(4.9) 



We display in Figure[8]the sensitivity functions with the superposed optimal grid obtained as in section 3.2.2 



Note that the maxima of the sensitivity functions are very close to the optimal grid points in the full 
measurements case. 

4.2.3 The preconditioned Gauss-Newton iteration 

Since the reconstruction mapping Q„ gives good reconstructions when properly interpolated, we can think of 
it as an approximate inverse of the forward map M„(Ao.) and use it as a non-Zmear preconditioner. Instead 
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of minimizing the misfit in the data, we solve the optimization problem 



min||Q„(vec(M„(A^))) - Q„(vec(M„(A^J))||^. 

(T>0 



(4.10) 



Here cr* is the conductivity that we would like to recover. For simplicity the minimization (4.f 1 is formulated 



with noiseless data and no regularization. We refer to |:17' for a study of the effect of noise and regularization 



on the minimization (4.101 



The positivity constraints in (4.10) can be dealt with by solving for the log-conductivity k — In(cr) 



instead of the conductivity cr. With this change of variable, the residual in (4.10) can be minimized with 



the standard Gauss-Newton iteration, which we write in terms of the sensitivity functions (4.5) evaluated at 
exp K^^': 



rU) 



^0+1) ^ ^0) _ ('diag(l/7(°') D^-f diag(e 



,U}' 



t 



(^expK^'j ) [Q„(vec(M„(A^^p^(,)))) - Q„(vec (M„(A^J))] . 

(4.11) 

The superscript f denotes the Moore-Penrose pseudoinverse and the division is understood componentwise. 
We take as initial guess the log-conductivity k'-'^' = Incr'^-', where a'^^'' is given by the linear interpolation 



of Q„(vec (M„(A(j^)) on the optimal grid (i.e. the reconstruction from section 4.2.1). Having such a good 
initial guess helps with the convergence of the Gauss-Newton iteration. Our numerical experiments indicate 



that the residual in (4.10) is mostly reduced in the first iteration [T3]. Subsequent iterations do not change 



significantly the reconstructions and result in negligible reductions of the residual in (4.10). Thus, for 



all practical purposes, the preconditioned problem is linear. We have also observed in ^} E] that the 
conditioning of the linearized problem is significantly reduced by the preconditioner Q„ . 

Remark 6. The conductivity obtained after one step of the Gauss-Newton iteration is in the span of the 



sensitivity functions (4.5). The use of the sensitivity functions as an optimal parametrization of the unknown 
conductivity was studied in I1 11. Moreover, the same preconditioned Gauss-Newton idea was used in \37l for 
the inverse spectral problem of section \3.S\ 

We illustrate the improvement of the reconstructions with one Gauss-Newton step in Figure u\ (c) . If 
prior information about the conductivity is available, it can be added in the form of a regularization term 



in (4.10). An example using total variation regularization is given in |13) . 



5 Two dimensional media and partial boundary measurements 

In this section we consider the two dimensional KIT problem with partial boundary measurements. As 
mentioned in section [T] the boundary B is the union of the accessible subset Ba and the inaccessible subset 
Bj. The accessible boundary Ba may consist of one or multiple connected components. We assume that 
the inaccessible boundary is grounded, so the partial boundary measurements are a set of Cauchy data 
{u\ig , {an ■ VM)|g }, where u satisfies (1.1) and u|g = 0. The inverse problem is to determine a from 
these Cauchy data. 

Our inversion method described in the previous sections extends to the partial boundary measurement 
setup. But there is a significant difference concerning the definition of the optimal grids. The tensor product 
grids considered so far are essentially one dimensional, and they rely on the rotational invariance of the 
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problem for a^°' = 1. This invariance does not hold for the partial boundary measurements, so new ideas 
are needed to define the optimal grids. We present two approaches in sections [5. 1| and [5?2l The first one uses 
circular planar networks with the same topology as before, and mappings that take uniformly distributed 
points on B to points on the accessible boundary Ba- The second one uses networks with topologies designed 
specifically for the partial boundary measurement setups. The underlying two dimensional optimal grids are 
defined with sensitivity functions. 

5.1 Coordinate transformations for the partial data problem 

The idea of the approach described in this section is to map the partial data problem to one with full 
measurements at equidistant points, where we know from section [4] how to define the optimal grids. Since 
ri is a unit disk, we can do this with diffeomorphisms of the unit disk to itself. 



Let us denote such a diffeomorphism by F and its inverse F ^ by G. If the potential u satisfies (l.ll, 



then the transformed potential u{x) = u{F{x)) solves the same equation with conductivity a defined by 

(5.1) 



a{x) = ^'^y>^y^ (^'(^)) 



T 



|detG'(y)| 



V=F{x) 



where G" denotes the Jacobian of G. The conductivity a is the push forward of a by G, and it is denoted 
by G*cr. Note that if G'{y) {G'{y)) ^ I and det G'{y) ^ 0, then ct is a symmetric positive definite tensor. 
If its eigenvalues are distinct, then the push forward of an isotropic conductivity is anisotropic. 

The push forward g*A^ of the DtN map is written in terms of the restrictions of diffeomorphisms G and 
F to the boundary. We call these restrictions g = G|g and / = -Fig and write 

{{g.A,)u,m - (A^K o g))(T)|,^/(,) , 9 e [0, 2tt), (5.2) 

for Mg G H^'^{B). It is shown in J64] that the DtN map is invariant under the push forward in the following 
sense 

g^K„ = Kca- (5.3) 



Therefore, given (5.3) we can compute the push forward of the DtN map, solve the inverse problem with 



data 5*Ao. to obtain G*ct, and then map it back using the inverse of (5.2 1. This requires the full knowledge 
of the DtN map. However, if we use the discrete analogue of the above procedure, we can transform the 
discrete measurements of A^ on Ba to discrete measurements at equidistant points on B, from which we can 
estimate a as described in section |4l 

There is a major obstacle to this procedure: The EIT problem is uniquely solvable just for isotropic 
conductivities. Anisotropic conductivities are determined by the DtN map only up to a boundary-preserving 



diffeomorphism [6l]. Two distinct approaches to overcome this obstacle are described in sections 5.1.1 and 
|5.1.2| The first one uses conformal mappings F and G, which preserve the isotropy of the conductivity, at 
the expense of rigid placement of the measurement points. The second approach uses extremal quasicon- 
formal mappings that minimize the artificial anisotropy of a introduced by the placement at our will of the 
measurement points in Ba- 
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Figure 9: The optimal grid in the unit disk (left) and its image under the conformal mapping F (right). 
Primary grid lines are solid black, dual grid lines are dotted black. Boundary grid nodes: primary x, dual 
o. The accessible boundary segment Ba is shown in solid red. 



5.1.1 Conformal mappings 

The push forward Gi^a of an isotropic a is isotropic if G and F satisfy G" [{G')'^^ 
This means that the diffeomorphism is conformal and the push forward is simply 



/ and F' {{F'f 



G*cr = (7 o F. 



(5.4) 



Since all conformal mappings of the unit disk to itself belong to the family of Mobius transforms [48], F 
must be of the form 



F{z) 



1 



z e 



UI<1, 



e[0,27r), aGC, |a| < 1, 



(5.5) 



where we associate M? with the complex plane C Note that the group of transformations (5.5) is extremely 
rigid, its only degrees of freedom being the numerical parameters a and w. 

To use the full data discrete inversion procedure from section |4] we require that G maps the accessible 
boundary segment Ba — {e*'^ | r € [— /3,/3]} to the whole boundary with the exception of one segment 
between the equidistant measurement points 9k, k — (n + l)/2, (n + 3)/2 as shown in Figure [9J This 



determines completely the values of the parameters a and oj in ( 5.5 ) which in turn determine the mapping / 



on the boundary. Thus, we have no further control over the positioning of the measurement points r^ = /(^fc), 
k = 1, . . . ,n. 

As shown in Figure l9] the lack of control over t^ leads to a grid that is highly non-uniform in angle. In 
fact it is demonstrated in |15| that as n increases there is no asymptotic refinement of the grid away from 
the center of Ba, where the points accumulate. However, since the limit n — >■ oo is unattainable in practice 
due to the severe ill-conditioning of the problem, the grids obtained by conformal mapping can still be useful 
in practical inversion. We show reconstructions with these grids in section |5.3[ 
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5.1.2 Extremal quasiconformal mappings 

To overcome the issues with confornial mappings that arise due to the inherent rigidity of the group of 

conformal automorphisms of the unit disk, we use here quasiconformal mappings. A quasiconformal mapping 

F obeys a Beltrami equation in O 

dF dF 

^=Ai(^)^- IIa^IIoo < 1, (5.6) 

with a Beltrami coefficient fi{z) that measures how much F differs from a conformal mapping. If /i = 0, 



then (5.6) reduces to the Cauchy-Riemann equation and F is conformal. The magnitude of ji also provides 



a measure of the anisotropy n of the push forward of a by F. The definition of the anisotropy is 



.(^..,.)^^^^H^, (5.7) 

v/Ai(z)/A2(z) + l' 

where Ai(z), X2{z) are the largest and the smallest eigenvalues of F^,a respectively. The connection between 
fi and K is given by 

KiF^,a,z)^\^l{z)\, (5.8) 

and the maximum anisotropy is 

K{F,<j)^svipK{F,<j,z)^\\n\\^. (5.9) 

z 

Since the unknown conductivity is isotropic, we would like to minimize the amount of artificial anisotropy 
that we introduce into the reconstruction by using F . This can be done with extremal quasiconformal map- 
pings, which minimize ||/i||co under constraints that fix / = -F|g, thus allowing us to control the positioning 
of the measurement points r^ — f{9k), for k — 1, . . . ,n. 

For sufficiently regular boundary values / there exists a unique extremal quasiconformal mapping that 
is known to be of a Teichmiiller type [53] . Its Beltrami coefficient satisfies 



M^)-IImII=o^, (5.10) 

for some holomorphic function 4>{z) in fi. Similarly, we can define the Beltrami coefficient for G, using a 
holomorphic function ip. It is established in [62^ that F admits a decomposition 

i^ = ^"1 o Aa' o $, (5.11) 

where 



^z) = j ^/W)dz, -^{O = J VW)dC, (5.12) 

are conformal away from the zeros of </> and ip, and 

AK{x + iy) ^ Kx + iy (5.13) 
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Figure 10: Teichmiiller mapping decomposed into conformal mappings $ and ^, and an affine transform 
Ak- The poles of and ^ and their images under $ and ^ are -^, the zeros of 4> and V' and their images 
under $ and ^ are D. 





Figure 11: The optimal grid under the quasiconformal Teichmiiller mappings F with different K. Left: 
K = 0.8 (smaller anisotropy); right: K — 0.66 (higher anisotropy). Primary grid lines are solid black, dual 
grid lines are dotted black. Boundary grid nodes: primary x, dual o. The accessible boundary segment Ba 
is shown in solid red. 



is an afhne stretch, the only source of anisotropy in (5.11 ) 



K{F,a) = M, 



K~ 1 



K^\ 



(5.14) 



Since only the behavior of / at the measurement points 9k is of interest to us, it is possible to construct 
explicitly the mappings $ and \1/ [15]. They are Schwartz-Christoffel conformal mappings of the unit disk to 



polygons of special form, as shown in Figure 10 See [151 Section 3.4] for more details. 

We demonstrate the behavior of the optimal grids under the extremal quasiconformal mappings in Figure 
[Tl] We present the results for two different values of the afhne stretching constant K. As we increase the 
amount of anisotropy from K = 0.8 to K — 0.66, the distribution of the grid nodes becomes more uniform. 
The price to pay for this more uniform grid is an increased amount of artificial anisotropy, which may 
detriment the quality of the reconstruction, as shown in the numerical examples in section |5.3[ 
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Figure 12: Pyramidal networks r„ for n 
X and the interior nodes with o. 



6, 7. The boundary nodes Vj, j — 1, . . . ,n are indicated with 



5.2 Special network topologies for the partial data problem 

The Umitations of the construction of the optimal grids with coordinate transformations can be attributed 
to the fact that there is no non-singular mapping between the full boundary B and its proper subset Ba- 
Here we describe an alternative approach, that avoids these limitations by considering networks with dif- 
ferent topologies, constructed specifically for the partial measurement setups. The one-sided case, with the 



accessible boundary Ba consisting of one connected segment, is in section 5.2.1 The two sided case, with Ba 



the union of two disjoint segments, is in section 5.2.2 The optimal grids are constructed using the sensitivity 
analysis of the discrete and continuum problems, as explained in sections [4.2.2| and |5.2.3[ 



5.2.1 Pyramidal networks for the one-sided problem 

We consider here the case of Ba consisting of one connected segment of the boundary. The goal is to choose 
a topology of the resistor network based on the flow properties of the continuum partial data problem. 
Explicitly, we observe that since the potential excitation is supported on Ba, the resulting currents should 
not penetrate deep into fl, away from Ba- The currents are so small sufficiently far away from Ba that in the 
discrete (network) setting we can ask that there is no flow escaping the associated nodes. Therefore, these 
nodes are interior ones. A suitable choice of networks that satisfy such conditions was proposed in [TB]. We 
call them pyramidal and denote their graphs by r„, with n the number of boundary nodes. 

We illustrate two pyramidal graphs in Figure [121 ^'-'^ n = 6 and 7. Note that it is not necessary that n 
be odd for the pyramidal graphs r„ to be critical, as was the case in the previous sections. In what follows 
we refer to the edges of r„ as vertical or horizontal according to their orientation in Figure [T2[ Unlike 
the circular networks in which all the boundary nodes are in a sense adjacent, there is a gap between the 
boundary nodes vi and u„ of a pyramidal network. This gap is formed by the bottommost n ~ 2 interior 
nodes that enforce the condition of zero normal flux, the approximation of the lack of penetration of currents 
away from Ba- 

It is known from |231 116) that the pyramidal networks are critical and thus uniquely recoverable from the 
DtN map. Similar to the circular network case, pyramidal networks can be recovered using a layer peeling 
algorithm in a finite number of algebraic operations. We recall such an algorithm below, from [IB], in the 
case of even n = 2m. A similar procedure can also be used for odd n. 
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Algorithm 1. To determine the conductance 7 of the pyramidal network (r„,7) from the DtN map A'"-*, 
perform the following steps: 

(1) To compute the conductances of horizontal and vertical edges emanating from the boundary node Vp, 
for each p ^ 1, • • • , '2m, define the following sets: 

Z = {vi,...,Vp^i,Vp+i,.. .,Vm}, C = {Wm+2,---,t'2m}, 

H = {wi, . . . , Wp} and V = {vp, . . . , Vm+i}, in the case p < m. 

Z = {Vm+l: ■ ■ ■ ,Vp^i,Vp+i, . . . ,V2m}, C = {vi, . . . ,Vm-l}, 

H = {vp, . . . , V2rn} and V = {«„, . . . , Vp}, form + l<p< 2m. 

(2) Compute the conductance ^{Epji) of the horizontal edge emanating from Vp using 



,(") 



^Z.H I ^H, 



compute the conductance ^{Epy) of the vertical edge emanating from Vp using 

^{Ep,^) = [a^}, - A^;:l (Ag),)" A^,%^ ly, 
where ly and 1h are column vectors of all ones. 



(5.15) 



(5.16) 



(3) Once ^{Ep,h), ^{Ep,^,) are known, peel the outer layer from r„ to obtain the subgraph r„_2 with the set 
S = {wi, . . . , W2m-2} of boundary nodes. Assemble the blocks K55, K^g, Kg^, Kg^ of the Kirchhoff 
matrix of (r„,7), and compute the updated DtN map A^"^^' of the smaller network (r„_2,7), as 
follows 



A("-2) = -K' ., - K, 



ss 



^SB 



(P (A(") - Kee) P' 



PK 



BS- 



(5.17) 



Here P G M'" 2)xn j^^ a projection operator: PP = ln-2, andK.'^^ is apart 0/K55 that only includes 
the contributions from the edges connecting S to B. 

(4) If m — I terminate. Otherwise, decrease m by 1, update n — 2m and go back to step 1. 

Similar to the layer peeling method in [22], Algorithm IT] is based on the construction of special solutions. 
In steps 1 and 2 the special solutions are constructed implicitly, to enforce a unit potential drop on edges Ep^ 
and Ep^^ emanating from the boundary node Vp. Since the DtN map is known, so is the current at Vp, which 
equals to the conductance of an edge due to a unit potential drop on that edge. Once the conductances are 
determined for all the edges adjacent to the boundary, the layer of edges is peeled off and the DtN map of a 
smaller network r„_2 is computed in step 3. After m layers have been peeled off, the network is completely 
recovered. The algorithm is studied in detail in [IB], where it is also shown that all matrices that are inverted 



in (5.151, (5.16) and (5.171 are non-singular 



Remark 7. The DtN update formula [5.11) provides an interesting connection to the layered case. It can 



be viewed as a matrix generalization of the continued fraction representation (3.36). The difference between 



the two formulas is that [3.36) expresses the eigenvalues of the DtN map, while (5.11) gives an expression 
for the DtN map itself. 
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Figure 13: Two-sided network r„ for n — 10. Boundary nodes Vj, j — 1, . . . ,n are x, interior nodes are o. 

5.2.2 The two-sided problem 

We call the problem two-sided when the accessible boundary Ba consists of two disjoint segments of B. A 
suitable network topology for this setting was introduced in jTTj. We call these networks two-sided and 
denote their graphs by T„, with n the number of boundary nodes assumed even n = 2m. Half of the nodes 
are on one segment of the boundary and the other half on the other, as illustrated in Figure [13] Similar to 
the one-sided case, the two groups of m boundary nodes are separated by the outermost interior nodes, which 
model the lack of penetration of currents away from the accessible boundary segments. One can verify that 
the two-sided network is critical, and thus it can be uniquely recovered from the DtN map by the Algorithm 
[2] introduced in [T7]. 

When referring to either the horizontal or vertical edges of a two sided network, we use their orientation 
in Figure [TSJ 

Algorithm 2. To determine the conductance 7 of the two-sided network (T„,7) from the DtN map Ay, 
perform the following steps: 

(1) Peel the lower layer of horizontal resistors: 

For p = 771 + 2, ni + 3, . . . , 2m define the sets Z = {p + l,p -\- 2, . . . ,p + m — 1} and C — {p ~ 2,p — 
3, . . . ,p — m}. The conductance of the edge Epqji between Vp and Vg, q = p — 1 is given by 



liEp,g,h) = -Ap,g -\- Ap,c(Az,c) ^^z,q 



(5.18) 



Assemble a symmetric tridiagonal matrix A with off-diagonal entries —"/(Ep^p^i^h) and rows summing 
to zero. Update the lower right m-by-m block of the DtN map by subtracting A from it. 

(2) Let s^m~l. 

(3) Peel the top and bottom layers of vertical resistors: 

For p = 1,2,..., 2m define the sets L — {p — l,p — 2, . . . ,p — s} and R = {p + l,p + 2, . . . ,p + s}. If 
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p < m/2 for the top layer, or p > 3m/2 for the bottom layer, set Z = L, C = R. Otherwise let Z ^ R, 
C = L. The conductance of the vertical edge emanating from Vp is given by 

l{Ep,v) = Ap_p - Xp^c{^z,cY^J^z,p- (5.19) 

Let D — diag('y(Ep^y)) and update the DtN map 

A^ = -D-D{A^ + Dy^D. (5.20) 

(4) If s — I go to step (7). Otherwise decrease s by 2. 

(5) Peel the top and bottom layers of horizontal resistors: 

For p — 1, 2, ... , 2m define the sets L = {p—l,p — 2, . . . ,p— s} and _R={p + 2,p + 3, ...,p + s + l}. // 
p < m/2 for the top layer, or p < 3m /2 for the bottom layer, set Z ~ L, C = R, q = p + 1. Otherwise 
let Z = R, C = L, q = p — 1. The conductance of the edge connecting Vp and Vq is given by (5.18). 
Update the upper left and lower right blocks of the DtN map as in step (1). 

(6) If s — Q go to step (7), otherwise go to (3). 

(7) Determine the last layer of resistors. If m is odd the remaining vertical resistors are the diagonal 
entries of the DtN map. If m is even, the remaining resistors are horizontal. The leftmost of the 



remaining horizontal resistors 7(-Ei,2.h) is determined from (5.18) with p = 1, q = m + \, C ^ {1, 2}, 



Z = {to + 1, to + 2} and a change of sign. The rest are determined by 

l{Ep^P+i,h) = (Ap,ff - Ap,c(Az,c)"'A2,ff) 1, (5.21) 

where p = 2,3,...,to— 1, C = {p— l,p,p+l}, Z = {p+to,— 1,p+to,,p+to + 1}, H = {p+to,— 1,p+7tt,}, 
and 1. is a vector (1, 1)"^. 

Similar to Algorithm [l] Algorithm [2] is based on the construction of special solutions examined in P^I24| . 
These solutions are designed to localize the flow on the outermost edges, whose conductance we determine 



first. In particular, formulas (5.18) and (5.19) are known as the "boundary edge" and "boundary spike" 



formulas O Corollaries 3.15 and 3.16]. 

5.2.3 Sensitivity grids for pyramidal and two-sided networks 

The underlying grids of the pyramidal and two-sided networks are truly two dimensional, and they cannot 
be constructed explicitly as in section [3] by reducing the problem to a one dimensional one. We define 
the grids with the sensitivity function approach described in section [4. 2. 2[ The computed sensitivity grid 
points are presented in Figure [M} and we observe a few important properties. First, the neighboring points 
corresponding to the same type of resistors (vertical or horizontal) form rather regular virtual quadrilaterals. 
Second, the points corresponding to different types of resistors interlace in the sense of lying inside the 
virtual quadrilaterals formed by the neighboring points of the other type. Finally, while there is some 
refinement near the accessible boundary (more pronounced in the two-sided case), the grids remain quite 
uniform throughout the covered portion of the domain. 
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Figure 14: Sensitivity optimal grids in the unit disk for the pyramidal network r„ (left) and the two-sided 
network T„ (right) with n = 16. The accessible boundary segments Ba are solid red. Blue D correspond to 
vertical edges, red -^ correspond to horizontal edges, measurement points are black x. 



Note from Figure 13 that the graph T„ lacks the upside-down symmetry. Thus, it is possible to come 
up with two sets of optimal grid nodes, by fitting the measured DtN map M„(Ao.) once with a two-sided 
network and the second time with the network turned upside-down. This way the number of nodes in the 
grid is essentially doubled, thus doubling the resolution of the reconstruction. However, this approach can 
only improve resolution in the direction transversal to the depth, as shown in [T71 Section 2.5]. 



5.3 Numerical results 

We present in this section numerical reconstructions with partial boundary measurements. The reconstruc- 
tions with the four methods from sections |5.1.1[ |5.1.2[ |5.2.1| and |5.2.2| are compared row by row in Figure 
15 We use the same two test conductivities as in Figure IWa). Each row in Figure 15 corresponds to one 
method. For each test conductivity, we show first the piecewise linear interpolation of the entries returned 



by the reconstruction mapping Q„, on the optimal grids (first and third column in Figure 15 1. Since these 



grids do not cover the entire i7, we display the results only in the subset of i7 populated by the grid points. 



We also show the reconstructions after one-step of the Gauss-Newton iteration (4.11) (second and fourth 



columns in Figure 15) 



As expected, the reconstructions with the conformal mapping grids are the worst. The highly non- 
uniform conformal mapping grids cannot capture the details of the conductivities away from the middle of 
the accessible boundary. The reconstructions with quasiconformal grids perform much better, capturing the 
details of the conductivities much more uniformly throughout the domain. Although the piecewise linear 
reconstructions Q„ have slight distortions in the geometry, these distortions are later removed by the first 
step of the Gauss-Newton iteration. The piecewise linear reconstructions with pyramidal and two-sided 
networks avoid the geometrical distortions of the quasiconformal case, but they are also improved after one 
step of the Gauss-Newton iteration. 

Note that while the Gauss-Newton step improves the geometry of the reconstructions, it also introduces 
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1.2 1.4 1.6 1.8 0.5 



Figure 15: Reconstructions with partial data. Same conductivities are used as in figure [71 Two leftmost 
columns: smooth conductivity. Two rightmost columns: piecewise constant chest phantom. Columns 1 
and 3: piecewise linear reconstructions. Columns 2 and 4: reconstructions after one step of Gauss-Newton 
iteration (4.11 ). Rows from top to bottom: conformal mapping, quasiconformal mapping, pyramidal network, 

Centers of supports of measurement (electrode) 



two-sided network 
functions Xq £^re x 



Accessible boundary Ba is solid red. 



41 



some spurious oscillations. This is more pronounced for the piecewise constant conductivity phantom (fourth 



column in Figure 15 1. To overcome this problem one may consider regularizing the Gauss-Newton iteration 



(4.11) by adding a penalty term of some sort. For example, for the piecewise constant phantom, we could 



penalize the total variation of the reconstruction, as was done in |13j . 

6 Summary 

We presented a discrete approach to the numerical solution of the inverse problem of electrical impedance 
tomography (FIT) in two dimensions. Due to the severe ill-posedness of the problem, it is desirable to 
parametrize the unknown conductivity cr(jK.) with as few parameters as possible, while still capturing the 
best attainable resolution of the reconstruction. To obtain such a parametrization, we used a discrete, model 
reduction formulation of the problem. The discrete models are resistor networks with special graphs. 

We described in detail the solvability of the model reduction problem. First, we showed that boundary 
measurements of the continuum Dirichlet to Neumann (DtN) map Acr for the unknown cr(x) define matrices 
that belong to the set of discrete DtN maps for resistor networks. Second, we described the types of network 
graphs appropriate for different measurement setups. By appropriate we mean those graphs that ensure 
unique recoverability of the network from its DtN map. Third, we showed how to determine the networks. 

We established that the key ingredient in the connection between the discrete model reduction problem 
(inverse problem for the network) and the continuum FIT problem is the optimal grid. The name optimal 
refers to the fact that finite volumes discretizations on these grids give spectrally accurate approximations 
of the DtN map, the data in FIT. We defined reconstructions of the conductivity using the optimal grids, 
and studied them in detail in three cases: (1) The case of layered media and full boundary measurements, 
where the problem can be reduced to one dimension via Fourier transforms. (2) The case of two dimensional 
media with measurement access to the entire boundary. (3) The case of two dimensional media with access 
to a subset of the boundary. 

We presented the available theory behind our inversion approach and illustrated its performance with 
numerical simulations. 
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A The quadrature formulas 



To understand definitions (3.8 1, recall Figure [I] Take for example the dual edge 



2 5j 2 ^ 2 '-^ 2 ' 



where Pj_i ,_i = fi(cos0j,sin0j). We have from (2.5|, and the change of variables to z{r) that 
o'(x)n(x) • Vw(x)(is(x) = 



rj_icr(ri_i) ^^ d9 '■ 



du{ri^i,dj) he (t/j-ij - U^,j) 



dr 



dz 



z{ri) - z(ri_i) 



which gives the first equation in (3.8). Similarly, the flux across 



is given by 



'iJ+h 



CT(x)n(x) • Vu(x)(is(x) = 



- (-Pi-ij+i'^^+ij+i)' 



2.J+2' ■'+2'JT2' 



. r 86 do 

Ut.j+i - U{i,j) 



(z(f,)-z(f,_i)) 



which gives the second equation in (3.8). 



B Continued fraction representation 

Let us begin with the system of equations satisfied by the potential Vj , which we rewrite as 



bj = bj+i + aj+iXVj+i, j = 0,1 



, . . . -u, 



bo = '^^ 

Vi+i = 0, 



where we let 



Vj ^Vj+i+a^bj. 



Combining the first equation in (2.2) with (2.2), we obtain the recursive relation 



1 



V, 



1 



J -1,2 



, . . . , ^, 



Sj+iA + 



bj+i 



V, 



j+i 



(2.1) 



(2.2) 



(2.3) 
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which wc iterate for j decreasing from j = £ — 1 to 1, and starting with 



be _ 1 

Ve ae 



(2.4) 



The latter fohows from the first equation in (2.1 1 evaluated at j = £, and boundary condition Ve+i — 0. We 
obtain that 



Ft(A) = l/i/$, = -i = 



Vi Vi 



1 



bo bi+aiWi jj^ 

ai\+ — 



has the continued fraction representation (3.361 



(2.5) 



C Derivation of results (3.47-3.48) 



^-ike 



To derive equation (3.47) let us begin with the Fourier series of the electrode functions 

fcGZ feGZ 

where the bar denotes complex conjugate and the coefficients are 



2ir 



C,{0) = ^j^ X,{6y''de = ^sinc 1^^ 



Then, we have 



Jo I. i,,cy Jo 



k.k'eZ 



-E' 



2tt 



,jfc(ep-e,)j/^2 



khf, 



PT^q- 



The diagonal entries are 



(^7)p.p--E(^7),, = -^E^^''''/(fc') 



9#P 



feez 



khf, 



Ee- 



ik6„ 



But 



Ee- 

q¥=p 



ikdf, 



Y^ -i^(g-l) _ -ikOp _ iTTk{l-l/n) sin(7rfc) ^-ikB^ _ „^_ JkB^ 



9=1 



sin(7rfc/r 



- e p = 7idfe.o - e p. 



(3.6) 



(3.7) 



(3.8) 



(3.9) 



(3.10) 



Since /(O) = 0, we obtain from ( [3J1|3.12[ ) that ([3J8]) holds for p = q, as well. This is the result ( |3.47 l 



Moreover, (3.48) follows from 



n 



9=1 



feiez 



sine 



fci/ie 



2 n 



Ee 

9=1 



z(/c— A;i)^g 



(3.11) 
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and the identity 



, 2,r(fc-fci) , 



^e.(fc-fci)«, = ^g.i^i^^F^(?-i) = „4^^^. 



(3.12) 



9=1 



9=1 



D Rational interpolation and Euclidean division 



Consider the case mi/2 = 1, where F{X) — 1/_F^^(A) follows from (3.36). We rename the coefficients as 



and let A = x^ to obtain 



K2J-1 = (ij, K2j = aj, j = 1,. . .£, 
F{x^) 



KlX 



(4.13) 



(4.14) 



K2X 



1 



K2i-lX + 



K2lX 



To determine Hj, for j = 1,. . . ,2i, we write first (4.14) as the ratio of two polynonfials of x, P2e{x) and 
Q2t-i{x) of degrees 21 and 2^—1 respectively, and seek their coefiicients Cj, 



Fjx"^) ^ P2i{x) ^ C2tx^^ + C2(£-i)a;^(^-^^) + . . . + C2X^ + cp 

X Qlf.^\{x) C2l-lX^^-'^ + C2e-3X^'^-^ + . . . + ClX 



(4.15) 



We normalize the ratio by setting cq = —1. 

Now suppose that we have measurements of i^ at A^ = x^, for fc = 1, . . . , 2£, and introduce the notation 



Xk 



^D, 



We obtain from (4.15) the following linear system of equations for the coefficients 



P2iixk) ~ DkQ2t~iixk) =0, k = l,...,2e, 



(4.16) 



(4.17) 



or in matrix form 



/ -Dixi 


x\ 


-Dixl .. 


. -Dixf-i 


^2e 

Xl 


-D2X2 


xl 


-D2XI . . 


. -D2xl'-' 


^21 
X2 


\ -D21X21 


^2 

■'^21 


-D2exlg . . 


. -D2ixff' 


^21 



/ Cl \ 



C2 



(4.18) 



V C.21 J 



with right hand side a vector of all ones. The coefficients are obtained by inverting the Vandermonde-like 



matrix in (4.18). In the special case of the rational interpolation (3.43), it is precisely a Vandermonde matrix 



Since the condition number of such matrices grows exponentially with their size 33 , the determination of 
{cj}j=i,...,2t is an ill-posed problem, as stated in Remarkfl] 

Once we have determined the polynomials P2i{x) and Q2e.-i{x), we can obtain {kj}j=i^...21 by Euclidean 



45 



polynomial division. Explicitly, let us introduce a new polynomial P2i-2{x) = C21-2X 



2(?-2 



. Co , so that 



H2X 



Q2e-i(x) 

P2i-2{x) 



KlX 



P2t-2{x) P2i{x) 



i2l-l 



{x) Q2l-l{x)' 



K3X- 



K2j'-ia; + 



K2ex 



Equating powers of x we get 



Kl = 



C2I 
C2£-l ■ 



and the coefficients of the polynomial P2t-2{x) are determined by 

C2j = C2j - K1C2J-I, j = 1,. . . ,(- 1, 

Co = Cq. 

Then, we proceed similarly to get K2, and introduce a new polynomial Q2i-3ix) so that 

1 P2£-2ix) , Q2i-3{x) Q2l-l{x) 



KsX- 



1 



Q2l-3{x) 



K2X 



P2l-2{x) P2i-2{x) 



KiX - 



1 



K2f-ia; + 



K,2eX 



Equating powers of x we get K2 = C21-1 / C21-2 and the polynomial Q2i-3(x) and so on. 



(4.19) 



(4.20) 



(4.21) 
(4.22) 



(4.23) 



E The Lanczos iteration 



Let us write the Jacobi matrix (3.31) as 



( -ai 61 ... 
bi -02 ^2 



\ 





/ 



V &,_i 

where —aj are the negative diagonal entries and hj the positive off-diagonal ones. Let also 

-A = diag(-^2 ...,-^2) 

be the diagonal matrix of the eigenvalues and 

Yj = diag(v/S^, . . . , v^)Yj 



(5.24) 



(5.25) 



(5.26) 
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the eigenvectors. They are orthonormal and the matrix Y = (Yi, . . . , Y^) is orthogonal 

YY^ = Y^Y = I. 
The spectral theorem gives that A = — YAY-^ or, equivalcntly, 



(5.27) 



AY = YA. 



(5.28) 



The Lanczos iteration [65|, I20j determines the entries aj and bj in A by taking equations (5.281 row by row. 
Let us denote the rows of Y by 

W, -ejY, j = !,...,£, (5.29) 



and observe from (5.281 that they are orthonormal 



W,W, - 5,- , 



(5.30) 



We get for j — \ that 






which determines Si, and we can set 



Wi = v/ar(ye 



The first row in equation (5.28) gives 



(5.31) 



(5.32) 



-aiWi+6iW2 = -WiA, 



(5.33) 



and using the orthogonality in (5.301, we obtain 



ai = WiAWf - Y. ^%^ ^1 = ll^iWi - WiAll 



and 



(5.34) 



W2 = 6r^ (aiWi - WiA) 



(5.35) 



The second row in equation (5.28) gives 



61 Wi - a2W2 + 62 W3 = -W2A, 



(5.36) 



and we can compute 02 and 62 as follows. 



a2=W2AW^, 62 = ||a2W2- W2A-61W1I 



(5.37) 
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Moreover, 



W3 = 62"! (02 W2 - W2A - 61 Wi) , 



(5.38) 



and the equation continues to the next row. 

Once we have determined {aj}j=i^...,i and {&j}i....,£-i with the Lanczos iteration described above, we 



can compute {aj,aj}j^i^,,,/. We already have from (5.31) that 






(5.39) 



The remaining parameters are determined from the identities 



a, = ^ dj^i + (1 - Oj,i) 



aitti 



Qfj \cij aj-i 



b,= 



"jV^ 



"J+1 



(5.40) 



F Proofs of Lemma [3] and Corollary [T] 

To prove Lemma 3 let A^^^ be the tridiagonal matrix with entries defined by {a^'' , Sj' }j=i,...^f, like in 
(3.29). It is the discretization of the operator in (3.58) with a ~-> a®. Similarly, let A*^") be the matrix 
defined by {ay' ,a° }j=i^,,,^i, the discretization of the second derivative operator for conductivity a'-°K By 



the uniqueness of solution of the inverse spectral problem and ( 3.82|3.83 ), the matrices A^"?) and A^"^ are 
related by 



diag 



:;(9) 






^Oi) 



(o) 



A^-^Miag 



:;(°) 



^3l^^'"--'\ 



a 



a 



(o) 
(9) 



= A(°)-aI. 



They have eigenvectors Y^'' and Yy' respectively, related by 



(6.41) 



diag 



Y 



(?) _ 



diag 



(o) 



J = !,...,£, 



(6.42) 



and the matrix Y with columns (6.42) is orthogonal. Thus, we have the identity 



(yy^) = s® J2 ^® = "1°^ E ^^"^ = 1' 



3 = 1 



j=l 



(6.43) 



which gives Si = a]^° by (3.82) or, equivalently 



:;(9) 



af' = ^^l = a®(0). 



Ao) 



(6.44) 



Moreover, straightforward algebraic manipulations of the equations in (6.41) and definitions (3.84) give the 



finite difference equations (3.85). D 
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To prove Corollary [l] recall the definitions (3.84) and (3.87) to write 



(6.45) 



Here we used the convergence result in Theorem [2] and denote by o(l) a negligible residual in the limit 
£ — >■ oo. We have 

7^®(C)rfC = Ea(°)a®(C(°))- / ''^'a®(C)dC + o(l), (6.46) 



''3 + 1 



and therefore 



S"+1 s+1 



<C 



'■j + 1 



a®(C)dC-5]S(°)a®(c(°)) 
p=i 



+ o(l), C=l/min(7(«)(C). 



(6.47) 



(o) 



But the first term in the bound is just the error of the quadrature on the optimal grid, with nodes at Q and 
weights a°' — C^\ ~ Q° , and it converges to zero by the properties of the optimal grid stated in Lemma 2 
and the smoothness of cr(^'(C). Thus, we have shown that 



^J + l ^3 + 1 



uniformly in j . The proof for the primary nodes C is similar. D. 



— ?► 0, as ^ — > oo. 



(9) 



(6.48) 



G Perturbation analysis 



It is shown in jl4t Appendix B] that the skew-symmetric matrix B given in (3.93) has eigenvalues HSj and 
eigenvectors 

where 



(7.49) 



(3^1 (<5,), . . . , yiiSj)f = diag (Sf , . . . , S|) Y„ (5^i(5,), . . . , 5^,(5,)) ^ = diag (af , . . . , a|) Y„ (7.50) 



Yj = {Yij, . . . ,Yi_j) are the eigenvectors of matrix A for eigenvalues —6^ and Yj = (Yi^ 
the vector with entries 



J , • ■ • , Ygj 



IS 



Y, 



P,3 






(7.51) 



G.l Discrete Gel'fand— Levitan formulation 



It is difficult to carry a precise perturbation analysis of the recursive Lanczos iteration that gives B from 

the spectral data. We use instead the following discrete Gel'fand-Levitan formulation due to Natterer [58] . 

Consider the "reference" matrix B'', for an arbitrary, but fixed r G [0, 1], and define the lower triangular, 
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transmutation matrix G, satisfying 



EGB = EB' G, ef G = ef 



(7.52) 



where E = I — e2£e^^. Clearly, if B = B*", then G = G'' = I, the identity. In general G is lower triangular 
and it is uniquely defined as shown with a Lanczos iteration argument in [141 Section 6.2]. 
Next, consider the initial value problem 



EB^(A) == iAE0(A), ef 0(A) = 1, 



(7.53) 



which has a unique solution 0(A) £ C^^, as shown in [TU Section 6.2]. When A ~ ±(5j, one of the eigenvalues 
of B, we have 

/ 2 



ct>i±s,) = y§-^yi±s,) 



"l^J 



■y{±6, 



and (7.53) holds even for E replaced by the identity matrix. The analogue of (7.531 for B*" is 

EB''0''(A) = iAE0''(A), ef ^'■(A) = 1, 



and, using (7.52) and the lower triangular structure of G, we obtain 



<P'-{±5,)^Gcj,{±S,), l<j<t 



Equivalently, in matrix form (7.56) and (7.54) give 



^"^ = G* = Gys, 



(7.54) 



(7.55) 



(7.56) 



(7.57) 



where $ is the matrix with columns (7.54), y is the orthogonal matrix of eigenvectors of B with columns 



( |7.49[ ), and S is the diagonal scaling matrix 

2 



o / ^ J- /'^-l/S ^-1/2 ,--1/2 /•- 1/2 

S = W— diagUi .^1 ,---,^i ,£.1 
\ ai \ 



Then, letting 

and using the orthogonality of y we get 



F = *'^S^ = Gy 



FF = GG^, 



where the bar denotes complex conjugate. Moreover, equation (7.52) gives 



EB'^F = EB'^Gy = EGB3^ = iEG^D = iEFD, 



(7.58) 

(7.59) 
(7.60) 

(7.61) 



where iG = zdiag ((5i, — (5i, . . . , (5^, —5t) is the matrix of the eigenvalues of B. 

The discrete Gel'fand-Levitan's inversion method proceeds as follows: Start with a known reference 
matrix B*", for some r € [0, 1]. The usual choice is B° = B'^°\ the matrix corresponding to the constant 
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coefRcient 17^°^ = 1. Determine $'' from (7.551, with a Lanczos iteration as explained in [TU Section 6.2]. 
Then, F = $''S^^ is determined by the spectral data (5'' and ^!^, for \ < j < L The matrix G is obtained 



from (7.60) by a Cholesky factorization, and B follows by solving (7.52), using a Lanczos iteration 



G.2 Perturbation estimate 

Consider the perturbations d5j = ASjdr and d£,j — A^jdr of the spectral data of reference matrix B*". We 
denote the perturbed quantities with a tilde as in 



D = D'' + rfD, s = s'^ + ds, y = y'' + dy, F = y + dF, 



(7.62) 



with D, S, y and F defined above. Note that F'^ = y , because C^ = I. Substituting (7.62| in (7.61) and 



using (7.55), we get 



EB'^dF = iEydB + ■^EdFD^ 



(7.63) 



T T T T 

Now multiply by y^ on the right and use that D''y — —iy^ B'' to obtain that dW = dFy satisfies 



EB'^dw - E dw B'' = iE ydB y 



with initial condition 



eJdW = eJdF y^ 



dJMl^dJ^'^' 



,dJ^,dJ^ r" 



(7.64) 



(7.65) 



Similarly, we get from (7.52) and C = I that 

EdB + 'EdGB'' ^FB'dG, ef dG = 0. 



Furthermore, equation (7.60) and F"" — y give 



dFy +ydF = dW + dW =dG + dG^. 



(7.66) 



(7.67) 



Equations (7.64 (7.67) allow us to estimate dfij/fi'^,. Indeed, consider the j,i + 1 component in (7.66) 



and use (7.67) and the structure of G, dG and B*" to get 

= dGj+i,j+i - dGj^j = dWj+i^j+i - dWjj, j = 1, . . . , 2^ - 1. 



df3, 



K 



(7.68) 



The right hand side is given by the components of dW satisfying ( 7.64||7.65 ) and calculated explicitly in [TH 
Appendix C] in terms of the eigenvalues and eigenvectors of B''. Then, the estimate 



2i-l 

E 



dp, 



/3!; 



<Cidr 



(7.69) 



which is equivalent to (3.98) follows after some calculation given in [TU Section 6.3], using the assumptions 

■J - 4°^ - rAS, and fj - ^f 



(3.74) on the asymptotic behavior of ASj and AS,j, i.e., of S^ — 5° = rA5j and f!^ — C° — fA(,j. 
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